Keywords: Acoustique sous-marine, inversion géoacoustique, méthode des sources images, acoustique des milieux multicouches, opérateur de Teager-Kaiser, représentations temps-fréquence underwater acoustics, geoacoustic inversion, image source method, acoustics of layered media, Teager-Kaiser energy operator, time-frequency representations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Ce travail de thèse propose d'analyser les signaux issus d'une source omnidirectionnelle sphérique réfléchis par un milieu sédimentaire stratifié et enregistré par une antenne d'hydrophones, en vue de caractériser quantitativement les sédiments marins aux moyennes fréquences, i.e. comprises entre 1 et 10 kHz.

La recherche développée dans ce manuscrit propose une méthodologie facilitant la recherche des paramètres géoacoustiques du milieu avec la méthode des sources images, ainsi qu'un ensemble de solutions techniques appropriées afin d'améliorer cette méthode d'inversion récemment développée.

La méthode des sources images repose sur une modélisation physique de la réflexion des ondes émises par une source sur un milieu stratifié sous l'approximation de Born. Par conséquent, la réflexion de l'onde sur le milieu stratifié peut être représentée par une collection de sources images, symétriques de la source réelle par rapport aux interfaces, dont les positions spatiales sont liées à la vitesse des ondes acoustiques et aux épaisseurs des couches.

L'étude se décline en deux volets : traitements des signaux et inversion des paramètres géoacoustiques.

La première partie du travail est focalisée sur le développement de la méthode des sources images. La méthode originelle se basait sur la construction de cartes de migration et de semblance de signaux pour déterminer les paramètres d'entrée de l'algorithme d'inversion qui sont des temps de trajet et des angles d'arrivée. Afin d'éviter cette étape, nous détectons les temps d'arrivée avec l'opérateur d'énergie de Teager-Kaiser (TKEO) et nous trouvons les angles par une méthode de triangulation. Le modèle d'inversion a été ensuite intégré en prenant en compte la possibilité de déformation de l'antenne. Cette partie se termine par une nouvelle approche qui combine TKEO et des méthodes temps fréquence afin d'avoir une bonne détection du temps d'arrivée dans le cas de signaux fortement bruités. Sur le plan du modèle et de l'inversion géoacoustique, nous proposons tout d'abord une description précise du modèle direct en introduisant le concept de sources images virtuelles. Cette étape permet de mieux comprendre l'approche développée. Ensuite, nous proposons une extension de la méthode des sources image pour l'inversion de paramètres géoacoustiques supplémentaires : la densité, l'atténuation et la vitesse des ondes de cisaillement. Cette extension est basée sur les résultats de l'inversion originelle (estimation iv du nombre de strates, de leur épaisseur, et de la vitesse des ondes de compression) ainsi que sur l'utilisation de l'amplitude des signaux réfléchis.

Ces améliorations et extensions de la méthode des sources images sont illustrées par leur application sur des signaux synthétiques et des signaux réels issus d'expérimentations en cuve et à la mer. Les résultats obtenus sont très satisfaisants, tant au niveau des performances de calcul que de la qualité des estimations fournies.

The image source method is based on a physical modeling of the wave reflection emitted from a source by a stratified medium under the Born approximation. As a result, the reflection of the wave on the layered medium can be represented by a set of image sources, symmetrical to the real source with respect to the interfaces, whose spatial positions are related to the sound speeds and the thicknesses of the layers.

The study consists of two parts : signal processing and inversion of geoacoustic parameters.

The first part of the work is focused on the development of the image source method. The original method was based on migration and semblance maps of the recorded signals to determine the input parameters of the inversion algorithm which are travel times and arrival angles. To avoid this step, we propose to determine the travel times with the Teager-Kaiser energy operator (TKEO) and the arrival angles are estimate with a triangulation approach. The inversion model is then integrated, taking into account the possible deformation of the antenna. This part concludes with a new approach that combines TKEO and time-frequency representations in order to have a good estimation of the travel times in the case of noisy signals.

For the modeling and geoacoustic inversion part, we propose first an accurate description of the forward model by introducing the concept of virtual image sources. This idea provides a deeper understanding of the developed approach. Then, we propose an extension of the image sources method to the estimation of supplementary geoacoustic parameters : the density, the absorption coefficient, and the shear wave sound speed. This extension is based on the results of the original inversion (estimation of the number of layers, their thicknesses, and the pressure sound speeds) and on the use of the amplitudes of the reflected signals.

Table des matières

List of figures

3.11

Uncertainties of sound speed (m/s) in each layer in the case of 5,10, 15 and 20 hydrophones for configuration 1 and seabed 1. . . . . . . . . . . . .

3.12

The effect of axial error detection on the image source position. . . . . . .

3.13

The effect of axial error detection on the sound speed profile . . . . . . . .

3.14

The effect of radial error detection on the image source position . . . . . .

3.15

The effect of radial error detection on the sound speed profile . . . . . . .

3.16

Error on sound speed in configuration 1 and seabed 1, in the case of axial and radial position errors in the third layer, with 5,10, 15 and 20 hydrophones

5.4

Representation of the distances between the source and the hydrophones and the distances between the first image source and the hydrophones . . 5.5 Shape of the antenna corrected with the triangulation and optimisation error approaches. Blue: the horizontal antenna, Black +: antenna shape found by the optimisation error approach, Green : antenna shape found by triangulation approach. Red: theoretical position. . . . . . . . . . . . . 5.6 Shape of the antenna corrected with the triangulation and optimized methods in the case of a perturbation in the order of three samples to each travel time detected. Blue: the horizontal antenna, Black +: antenna shape found by optimized approach, Green : antenna shape found by triangulation approach. Red: theoretical position. . . . . . . . . . . . . . . . . . .

Correction

Introduction

The underwater world covers 70% of the Earth surface and collecting information on this environment is essential to understand its huge natural resources. This makes the study of the ocean the center of economic and scientific issues. In fact, the knowledge of the seabed has become essential in many fields. For energetic economic interests, oil and hydrocarbons exploitation is highly related to the seafloor nature. From a military point of view, the seafloor is either the low boundary of the water column and thus determines the detection range of sonar systems, or the place where the mines are located and it is mandatory to have a good seafloor knowledge to achieve reliable mine detection. Geologists make an intensive use of acoustic data to discover the history of the Earth. In civil engineering, establishing structures such as offshore platforms, wind turbines, water turbines or settling cables depends on the seabed nature.

The multiplicity of the methods used by sedimentologists to achieve the characterization of the seabed is the translation of its complexity. To accomplish this necessary seabed knowledge, many devices have been deployed. They can be divided in two families. The first one consists of the instruments that make direct measurements in the seafloor or that grab sediment samples in order to make measurements in a laboratory. If the core drilling techniques or direct measurements provide valuable information about the seabed nature, they are still very expensive. Moreover, they give only information locally, at a precise location. It is sometimes tricky to extend the collected information to larger areas.

The second family is composed of systems that collect information remotely. They are based on different physics principles such as photography, magnetism, optics, or acoustics. They overcome the two difficulties pointed out previously (cost and localization) but they have of course their own drawbacks. Among these remote sensing techniques, the ones that use acoustics are the most developed because the sea water is a very good conductor of acoustic energy : in the ocean, the sound is able to propagate to hundreds or thousands of kilometers. Seabed characterization by acoustics can be achieved by different instruments or techniques. For example, side-scan sonars produce seafloor images with very high resolution from which the nature of the sediment interface can be inferred.

The geoacoustic inversion methods aim at estimating the properties of the seafloor by studying the impact of this environment on acoustic propagation. These methods might be active, using a controlled emission with a known acoustic source, or passive, using opportunity sources such as boats or marine mammals. The common idea of these methods is that, knowing the forward model of acoustic interaction with the seafloor, it might be possible to estimate seafloor parameters by analyzing the acoustic field reflected or scattered. Depending on the frequency and the acquisition geometry, various penetration depths and resolution are obtained. Similarly, different methods will provide estimation of different parameters. Among all these parameters, the most common are probably the number of layers that compose the layered seafloor, their thicknesses, and the compression wave sound speeds.

In this context, the work presented in this thesis is devoted to the improvement and extension of a recently developed geoacoustic inversion method, the image source method (ISM). ISM provides the local sound speed profile (SSP) from the recording by a hydrophones array of acoustic signals emitted by a point source and reflected by the seafloor. Working in the mid-frequency band (roughly between 1 and 10 kHz), it combines a good penetration depth, tens of meters depending on the seabed nature, and a good resolution, between 50 cm and 1 m. Due to its local properties, ISM is able to invert range-dependent environments, providing that the system is mounted on moving vehicles, for instance, on a AUV (Autonomous Underwater Vehicle) towing a source and an array.

The fundamental idea of the method is based on the original description of the reflection of spherical acoustic waves by stratified media by image sources whose positions are directly related to the SSP. The detection and localization of these sources gives the input data to estimate the SSP of the sedimentary structure. This forward model reduces drastically the number of inputs of the inverse problem and, consequently, the ISM has a low computational cost.

ISM can be described as a two-step algorithm : first, from the recorded signals, the image sources are detected and localized. Second, from the positions of these image sources, the thicknesses and sound speeds of each layer of the seabed are estimated by a recursive calculus simply based on the Snell-Descartes laws of refraction.

The major improvement of ISM done in this thesis concerns the detection and localization of image sources from the recorded signals. To do so, we propose to use the Teager-Kaiser energy operator. This operator gives a local estimation of the energy needed to produce a signal. In this work, we use its property of good spike detection. Therefore, the time delays between the image sources and the hydrophones can be estimated accurately and, after a triangulation process, the image sources are localized and the inversion algorithm can be performed.

Initially, the ISM provides only the SSP of compression wave inside the multilayered seafloor. In this work, we propose two extensions. First, we have shown on synthetic data that it is possible to use ISM to detect the presence of shear waves in the seafloor and to estimate their sound speed. Second, by using the amplitudes of recorded signals as a supplementary input, we have shown theoretically and on both synthetic and real data that it is possible to have an estimation of the densities of the sediment layers. This manuscript is organized into six chapters and three appendixes.

Chapter 1

The first chapter describes the sedimentary environment, which is the core of our study, and the geoacoustic inversion methods. The sediments are first defined and described. Then, their main properties are presented. Depending on the type of studies, different parameters might be used to characterize the sediments. Here, we focus on the properties related to the geoacoustic problem : the stratification, the classification, the porosity, the sound speed, the density, and the attenuation.

The second part of this chapter presents the geoacoustic inversion methods from a general point of view. Rather than describing in details some particular methods, we preferred to give the methodology generally used in these approaches in order to enlighten the similarities and discrepancies with ISM.

Chapter 2

The second chapter presents the signals on which ISM is applied in this work. Two types of signals are used. The first ones are synthetic signals which are produced in order to evaluate the performances of the ISM. The signals are obtained through a numerical evaluation of the Sommerfeld integral which enables to evaluate the reflection of a spherical wave by a complex seafloor. Two different geometric configurations and four different seafloors are used to produce these simulated signals.

The second type of signals are those which come from real experiments. We performed tank experiments in order to have signals obtained in a controlled environment and thus, to be able to compare the inversion results with known ground truth. The ISM is also tested on at-sea data coming from a campaign that took place in the Mediterranean Sea in 1998.

Chapter 3

In this chapter, the principle, algorithm, and sensitivity of ISM are studied. First, the representation with image sources of a reflected wave by a layered seafloor is detailed. A special focus is made on the distinction between "real" and "virtual" image sources which is one of the keys of the ISM. Then, assuming that the first step of ISM (detection and localization of the image sources) is already made, which is described in chapters 4 and 5, the estimation of SSP through Snell-Descartes laws is described.

In the last section of this chapter, the limits and sensitivity of ISM are described. On one hand, the study of the model error due to the equivalent virtual image sources with different type of seabed and different length of antenna is presented. On the other hand, the effect of radial and axial errors in detection of the equivalent virtual image sources are also studied.

Chapter 4

The work presented in the fourth chapter deals with the first step of the ISM which is the travel time detection. For this, in the first part of this chapter, the property of TKEO is exploited. This approach avoids the computation of semblance and migration maps and reduces the number of false alarms compared to the original algorithm. This study is completed by an improvement of the resolution and a strategy to detect and Introduction correct the presence or the absence of spikes in the detection matrix.

In the second part of this chapter, for more robustness to noise, the detection of peaks is done in time-scale or in time-frequency domains. The TKEO is applied to the coefficients at each scale or in each line and each column of spectrogram presentation. Then, the travel times detected by each approaches are compared to the theoretical travel times.

Chapter 5

The fifth chapter presents the triangulation algorithm to localize the equivalent virtual image sources. A study of the sensitivity of this algorithm is also presented. Duo to its high sensitivity to the array geometry, two corrections methods of the localization algorithm with advantage of fast calculations and easy implementation are presented.

Finally, the ISM is applied on the signals presented in chapter 2. The results are satisfactory. From a computational point of view, the method is faster than the original one based on migration and semblance maps. From an inversion point of view, the geoacoustic parameters are well estimated.

Chapter 6

In this chapter, for a more complete characterization of seabed, we have extended the ISM to the estimation of the sound speed of shear wave, the densities, and the absorption coefficients in sedimentary layers. First, the presentation of the new approach used to detect and estimate the sound speed of shear waves is detailed. Second, two different approaches are used to estimate the densities and the absorption coefficient of the stratified medium by using the parameters calculated by the ISM (the sound speed of pressure waves, thickness of layers and arrival angle) and the amplitudes of the reflected signals. These two methods, based on iterative estimations of geoacoustic parameters, are tested on the data presented in chapter 2. Finally, this approach is tested on the tank data.

Chapter 1

Sediments and Acoustics

Introduction

In oceanic engineering, identifying the nature of the seabed is of prime importance. Due to the difficulty and the high cost of sampling sediments directly from the seabed, the marine geology needs tools to perform a rapid analysis of the structure of a sedimentary area. Thus, many acoustic methods have been developed to identify the seabed through the properties of the acoustic wave which propagates in sedimentary environments.

The acoustic methods used for identifying the nature of the seabed have many potential applications, namely, economic, energetic, military, and scientific. The differences in application areas change the acoustic methods used to identify the seabed, such as the frequency band, the application mode, the resolution and the type of emitted signal.

CHAPTER 1. SEDIMENTS AND ACOUSTICS

In order to better understand these acoustic methods for seabed identification, this chapter is divided in two parts. The objective of the first part is to define the marine sediments geologically in order to identify some geoacoustic features that will be useful to describe the propagation of the sound wave in this medium. The second part presents the different types of acoustics systems and inversion algorithms used for estimating the different geoacoustic parameters of the marine sediments at different scales according to the instrumentation used.

Marine sediments

Definition

Marine sedimentation includes different phenomena such as erosion, friction, collisions, chemical dissolution, transport by currents, accumulation and deposition of abiotic particles. This also includes living or dead organic matter which form the marine snow and go down to the bottom under the effect of gravity to gradually form deposit layers that formed in the past most of the sedimentary layers, drawing changing sedimentary features. The formation of sediments leads, in most cases, to their vertical stratification. The geometry of the studied environment varies slowly in the horizontal plan but quickly in the vertical plan.

This thesis is focused mainly on the superficial seabed of the continental plateau (also called continental shelf). The continental shelf constitutes 20% of the seabed that is about 14% of the Earth's surface [START_REF] Pouliquen | Identification des fonds marins superficiels a l'aide des signaux d'echo-sondeurs[END_REF]. As its name suggests, this area is part of the continent ; it ensures the transition between the land surface and the deep ocean bottom. The depth of the bottom on the plateau is relatively low ; its evolution from the coast to the open sea is slow to about 200 m. The end of the continental plateau called "slope" is marked by a sudden and strong descent to the deep are as (Figure 1.1).

Structure of marine sediments

The structure of marine sediments provides important information on sedimentary environments such as the flow direction, the climate (arid and humid areas) and the setting (eg., river or sea). These structures may also indicate the origin of rocks and 1.2. MARINE SEDIMENTS 9 their destination (Figure 1.2). But there is no overall trend on the sequence of deposits, stipulating that the materials become less rude to landfill. However, the sedimentary rocks are often deposited in discrete layers, which lead to a particularly important sedimentary structure called bedding. . The bedding layers can vary in thickness from several millimeters to several tens of meters. Generally, but not always, the bedding is originally a horizontal orientation. An inclined bedding indicates that the rock has been, in any way, deformed. For this, some generic laws on sedimentary sequences were deduced by studying the deposit process [START_REF] Reineck | Depositional sedimentary environements with reference to terrigenous clastics. 2nd edition revised[END_REF], the environmental conditions, the stratigraphy [START_REF] Nichols | Sedimentology and stratigraphy[END_REF] and the observation of soil core samples. These laws are filed in compliance with the principles set out by Nicholas Steno [Nic99] :

-The principle of primary horizontality : the sedimentary layers are deposited horizontally. A sedimentary sequence which is not in horizontal position has undergone deformation subsequent to its deposition.

-The superposition principle : in the absence of structural changes, one layer is more recent than the one it covers and is older than the one that covers it. It is the oldest in geology.

-The principle of lateral continuity : the sedimentary layers are deposited laterally that is to say horizontally in all directions until the thickness becomes zero or they end on the banks of the sedimentary basin.

To summarize, the sediments are stratified and can be assimilated as a stack of homogeneous layers of different natures. This stack varies rapidly with depth, but may vary slowly in the horizontal direction. An example of a core made during the measurement campaign in the region of Bouma (Senegal) is presented in figure 1.3.

Geological property of marine sediments

There are many geotechnical parameters that characterize the sediment but in this thesis, we limit the study to the ones that affect most the acoustic problem. Therefore, we are only interested in the classification, the porosity and the evolution with depth. The acoustic parameters are presented in the next section.

Classification

The analysis of marine sedimentation has shown the heterogeneity and the complexity of marine sediments. Due to this diversity, the classification of the floor is a difficult task because it is reduced to the classification of a multitude of existing sediments into a limited number of categories. Despite these difficulties, two main classification strategies are usually used :

-mineralogical classification based on the chemical compounds of the grains : this classification aims at determining the chemical nature of the dominant components.

-particle size classification based on the grain size distribution and some related statistics such as the mean grain size, the skewness or the sorting.

A mineralogical classification is not suitable for the present study because it requires a very detailed analysis after sampling. The appropriate classification, used in this thesis, is the one based on the grain size, the particle size of which they are composed, given the strong bond of this parameter with the geoacoustic parameters that form the subject

GEOLOGICAL PROPERTY OF MARINE SEDIMENTS

of this study [START_REF] Jackson | High-Frequency seafloor acoustics[END_REF][START_REF] Tucker | Techniques in sedimentology[END_REF]. A commonly used classification, based on the size of the grains, is performed into four classes : gravel, sand, silt and clays. Note that clays and silts are often grouped under the term mud. The diameter of the grains is often given in logarithm base 2 by φ d = -ln(d) where d is the diameter in millimeters [START_REF] Jackson | High-Frequency seafloor acoustics[END_REF] (Table .1.1).

Gravel is composed of rock elements or inconsistent agglomerates of varied origin (eg., mollusk shells). These are particles with a bigger size than those of sand. Their size is between 4 and 20 mm (sometimes up to 50 mm). The interstices between the particles are filled with water or sand when the sediment is not unimodal [START_REF] Pouliquen | Identification des fonds marins superficiels a l'aide des signaux d'echo-sondeurs[END_REF].

The majority of sand formation has a detrital origin and is related to the rock erosion, which is reflected in its geological nature. It is a loose sediment whose elements have a size comprise between a few tens of micrometers to about 2 mm. The mineralogical nomenclature of sands is established according to the nature of the dominant elements : calcareous sands, quartz sands, micaceous sands and coral sands. The particle size nomenclature is quite variable according to the authors. For example, Sleath proposes a classification of sands as follows [START_REF] Sleath | Seabed mechanics[END_REF] :

-very coarse sand : 1 to 2 mm -coarse sand : 500 microns to 1 mm -medium-size sand : 250 to 500 microns In geology, silt is a sedimentary formation whose grains are intermediate in size between clay and sand that is to say between 4 and 50 micrometers (the exact limits can vary a bit according to laboratories). Silt characterizes the wind loess deposits, also called "silt plateaus" and is common in alluvial deposits. In the latter case, they have a major interest for the renewal of soil fertility and thus for agricultural work [START_REF] Sleath | Seabed mechanics[END_REF].

Historically, in geology and soil science, the term clay corresponds to all minerals having a size less than 2 microns. This particle size cut is inherited from petrographic studies by optical microscopy in the late nineteenth century. Then, the crystals having a size less than 2 micrometers were not recognizable and were classified under the designation clay. Today, the term clay differs depending on the areas of study. Thus, in geotechnics, where one is primarily interested in the mechanical behavior of soils, clay denotes the particle size materials below 4 micrometers (between 4 and 50 micrometers, it is called silt). In clay science, clay does not correspond to a particle size cut, but to minerals. The term is then used to describe the phyllosilicates and more particularly the clay minerals. The Wentworth classification is not enough to describe fully the sediment [START_REF] Garlan | Apports de la modélisation dans la sédimentologie marine récente[END_REF]. When the sediment cannot be described with a single type of grain, the Sheppard [START_REF] Jackson | High-Frequency seafloor acoustics[END_REF] ternary diagram allows a more detailed classification (Figure 1.4). At each end of the diagram there is a primary component of the sediment and the intermediate positions are mixtures of these constituents. Some example like sandy clay, silty sand or sandy slit are presented with their pourcent respectively. It should be noted that this presentation is only interested in soft sediment, it does not classify gravel, stones and pebbles [START_REF] Le | Comparaison de mesures de réflectivité à l'interface eau sediments à partir de données de sondeurs multifaisceau, de sédiments et de sismique Haute Résolution[END_REF].

Porosity

Porosity is a reference parameter for geologists and it is present in almost all geoacoustic models. The porosity is the set of empty spaces of a solid material ; these empty 1.4. GEOACOUSTIC PARAMETERS spaces are filled with a fluid (liquid or gas). Porosity is defined as the ratio between the volume of empty spaces and the total volume of a porous medium [START_REF] Cojean | Sédimentologie[END_REF] :

φ = V E V T , (1.1)
where V E is the volume of the fluid (water) and V T is the total volume of the sediment.

The porosity increases inversely with the size of the grains [START_REF] Hamilton | Sound velocity and related properties of marine sediments[END_REF]. In fine sediments like clay and sand, porosity is very high on the surface. Then, its value decreases rapidly over the first few meters, which proves that the seabed parameters evolve according to the applied pressure.

Gradients

The properties of marine deposits and geoacoustic parameters such as sound speed, density and attenuation change depending on the pressure applied with the landfill [START_REF] Hamilton | Variation of density and porosity with depth in deep-sea sediments[END_REF]. This is due to the porosity which decreases rapidly with depth for clays and silts (65% to 40% from 0 to 500 m) but varies slightly in a landfill up to 100 m for sandy bottoms, as they undergo very little reorganization and also because of their non-cohesive character. The change in porosity is due to compaction and cementation generated by the pressure which affects the sound speed and density which increases with immersion while mitigation decreases. The sound speed and density have gradients depending on the depth of the sediment layer. These gradients are about 1 m/s -1 per m for the sound speed and 1.4 10 -3 per m for the density. Generally, mechanical compaction and cementation are zero over the first hundred meters [START_REF] Cojean | Sédimentologie[END_REF]. However, in some cases, when the sound speed in the surface of the sediment is lower than that in the water or in the case of a graded bedding in a layer, this may induce the variations of acoustic parameters on a small thickness to the order of a few decimeters.

Geoacoustic parameters

The characterization of the seabed by acoustic requires deep knowledge of the acoustic behavior of sediments. A geoacoustic model of the soil must be build, i.e a model which link the geotechnical characteristics of a type of floor and its acoustic parameters. The sound speed (of compression P waves and shear S waves), the density and the attenuation coefficient are used by acoustician to identify the material according to its mechanical properties. But this remains insufficient because just two or three parameters can not really identify the sediment type. Therefore, these parameters are alternately examined in this report. The relationships between these parameters are complex, given the variety of these phenomena in the environment, such as gas bubbles or compaction and cementation of sediments and the sediment type. For example, the P and S waves in the lithified sediments are both comparable to that sound speed of water, whereas in the non lithified sediment, only P-waves are comparable to that of water.

In the literature, there are two approaches to get these settings : the first is based on a database from many works (see for example [START_REF] Hamilton | Variation of density and porosity with depth in deep-sea sediments[END_REF][START_REF] Hamilton | Compressional-wave attenuation in marine attenuation sediments[END_REF][START_REF] Hamilton | Reflection coefficients and bottom losses at normal incidence computed from pacific sediment properties[END_REF]) which describe the measured relationships between the geological and the acoustic parameters CHAPTER 1. SEDIMENTS AND ACOUSTICS of the widest possible variety of sediments. In general the approximation of fluid medium to describe these types of sediments is appropriate in the sense that the attenuation of S waves is high [START_REF] Hamilton | Variation of density and porosity with depth in deep-sea sediments[END_REF][START_REF] Jackson | High-Frequency seafloor acoustics[END_REF]. The advantage of this method is that it provides trends and ranges of possible values for these parameters and allows comparisons to be made between the different types of floor. On the other hand, the values given globally from point measurements are not accurate to treat individual cases.

The second method is based on a theoretical approach inspired from a physical and a mechanical description of the environment. Equations relate a model of wave propagation directly to geological and physical parameters. In this method, two theories are more used than others :

-the Biot-Stoll theory where the sediment is modeled by a porous medium composed of an isotropic elastic skeleton and a viscous fluid [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated pourous solid[END_REF],

-the Buckingham theory where the sediment is modeled by a granular unconsolidated medium saturated by a viscous fluid [START_REF] Buckingham | Theory of acoustic attenuation, dispersion, and pulse propagation in unconsolidated granular materials including marine sediments[END_REF].

Sound speed

In a fluid medium, the sound propagates only by compression. If the sound is not too loud, the compression and expansion of the fluid can be considered isentropic and its sound speed is the square root of the partial derivative of the acoustic pressure p by the density ρ s at constant entropy :

c f luid = dp dρ s (1.2)
The speed of sound in water is a function of the density ρ and the adiabatic compressibility coefficient χ :

c water = 1 ρχ (1.3)
As any porous or granular material, the sediments admit a certain elasticity, depending on their state of consolidation. Therefore, two types of waves can propagate : the compression waves P and the shear waves S. In the literature, the assumption of neglecting the S waves is highly used in the sense that the amplitude of these waves is negligible compared to the amplitude of P wave. Moreover their speed c s is often very low in sediment compared to the speed of compression wave c p and their attenuation coefficient is much more higher.

The sound speed is the main parameter to considered, for the simple reason that the greater part of geoacoustic parameters can be identified as a function of sound speed. For example, the speed of an acoustic wave in the sediment is related to the porosity and the density. It is possible to go back to the porosity and the particle size, from the sound speed values which are one of the most accessible parameters of sonar systems [START_REF] Robinson | Geophysical signal analysis[END_REF]. Hence, the sound speed values of the sediments have a clear interest in the knowledge of the soil. Hamilton and Bachman proposed a series of sound speed values in sediments and link them to different geotechnical parameters [START_REF] Hamilton | Sound velocity and related properties of marine sediments[END_REF]. Figure 1.5 shows the sound speed values depending on the porosity φ : the sound speed increases when porosity decreases. We also note that the sound travels more slowly in some sediments than in water. This particular property is observed in sediments where the porosity is very strong and the stiffness is very low [Gui99

Density

The density can be measured quite easily in the laboratory. It is the ratio of the mass density of the material over the mass density (in kg/m 3 ) of the water taken in the same conditions. Density is said to be apparent if the divisor volume is the total volume of the sample (wet sediment), and absolute if the divisor volume is the volume of the grains (dry sediment). The volume densities or the densities reported in this thesis are CHAPTER 1. SEDIMENTS AND ACOUSTICS those of the wet sediment. Figure 1.6 shows the density depending on the grain size, according to Hamilton and Bachman [START_REF] Hamilton | Sound velocity and related properties of marine sediments[END_REF]. Just like the sound speed, the density and its direct relation to porosity is a significant parameter of the sediment. As illustrated in figure 1.6, the density shows correlation with the grain size. There is a large dispersion in the values that can be explained by the shape of the grains, the percentage of fine particles and therefore a different porosity to the particles in the same size class. To assign a density to a given sediment, it is necessary to look in the literature for a value given for sediments having not only nearly a similar size, but also a porosity close to the considered sample. Like the porosity, the density, which is in direct relation with this parameter, can change as a result of the landfill. Few measurements have been made. The variation in density with the depth is analyzed for a configuration of the deeper bottom in [START_REF] Hamilton | Variation of density and porosity with depth in deep-sea sediments[END_REF]. Some values of the densities of different sedimentary layers are presented in table.1.2.

Acoustic impedance

Two types of impedance can be defined. The first one is the characteristic impedance of the propagation medium, defined by :

Z = ρc , (1.4)
where ρ is the density and c is the sound speed in the medium.

The second impedance is defined as the ratio of the acoustic pressure (p) over the 1.4. GEOACOUSTIC PARAMETERS particular velocity (v) :

Z = p v .
(1.5)

These two impedance are equal in the case of a progressive harmonic plane wave.

When an acoustic wave encounters an interface separating two media with different acoustic impedances, a part of the wave is transmitted into the other medium, while the other part is reflected (Figure 1.7). The acoustic impedance is a tool to study completely and quantitatively this phenomenon and to estimate the amount of acoustic energy transmitted and reflected.

Water (c 0 , ρ 0 ) Layer 1 (c 1 , ρ 1 ) θ (0) 1 θ (0) 0 θ (0) 0 R T Figure 1
.7 -Reflection and transmission of plane acoustic wave at a plane interface between two fluid media.

In the case of a harmonic progressive plane wave, the normal impedance of medium i is given by :

Z ni = p v z = ρ i c i cos(θ i ) , (1.6)
where v z is the projection on the axis normal to the interface (here, it is the (Oz) axis) of the particular velocity.

The boundary conditions on the interface for the pressure the normal velocity therefore consider that there is a continuity of the normal impedance on the interface :

Z n 1 = Z n 2 .
(1.7)

From this continuity condition, the coefficients of reflection and transmission are obtained :

R = Z n 1 -Z n 2 Z n 1 + Z n 2 , and, T = 2Z n 2 Z n 1 + Z n 2 (1.8)
The interest of the normal impedance is to study, for example, the propagation in the stratified medium. By using the signals reflected by a medium, we cannot always go back clearly to the sound speed and the density but to the reflection coefficient showing the impedance contrast between two layers [START_REF] Panda | Sediment classiffication based on impedance and attenuation estimation[END_REF][START_REF] Schimel | Contribution au dévéloppement de l'offshore pétrolier[END_REF]. So the acoustic impedance is an important parameter in acoustics to study the reflection and transmission of waves (Equation 1.8) in order to find the geoacoustic parameters.

CHAPTER 1. SEDIMENTS AND ACOUSTICS

Attenuation

Definition Attenuation is the decrease of the acoustic intensity or pressure level during the propagation of a sound wave in a medium between two points. It is caused by three independent phenomena : -the geometric divergence loss which is usually denoted by T L for transmission loss, quantifying how the acoustic energy is distributed in the space during the propagation of the wave. For the spherical wave, it is equal to T L = 20 log(r) and for the cylindrical wave it is equal to T L = 10 log(r) where r is the range. -the absorption which is the effect of the viscosity, the molecular relaxation, and the inter-grain friction within the sediment. -the interaction with obstacles (an interface or an object) in the water which generate scattering or reflection.

Absorption in sediments

Different units can be used to quantify the absorption of an acoustic wave in a medium, the main ones are presented in this paragraph. Note that the symbols used here are not standardized and might differ form one author to an other.

α in dB/m (lower level of intensity on 1 m), β in dB/λ (lower level of intensity on the wavelength λ), κ in dB/m/kHz (lower level of intensity on 1 m by kHz), -N e in Np/m (Neper by meter). The first three units quantify a decreased level for a distance while the latter applies to the amplitude of the pressure. If A 0 is the initial amplitude of the signal and A(r) is the received amplitude at distance r, N e is thus defined by the following relation : (-Ner) .

A(r) = A 0 e
(1.9)

The passage from one unit to an another is given in table The marine sediments are well known for their absorbent property [START_REF] Hamilton | Compressional-wave attenuation in marine attenuation sediments[END_REF] and the value of absorption discriminates the sediment type [LMR + 92]. The measurement principles of this parameter have been well studied for example in references [Ger04, TGS + 06, RWB + 02]. The absorption is connected to several parameters of the emitted signal such as the frequency as well as the nature of the ground grain sizes and porosity. Figure 1.8 shows the absorption of compression waves as a function of porosity and grain size. Figure 1.8 shows that the absorption increases with the grain size, up to about 50% for the porosity and nearly 4 for the mean grain size, beyond coming down. This border is the passage of silt to sand [START_REF] Le | Comparaison de mesures de réflectivité à l'interface eau sediments à partir de données de sondeurs multifaisceau, de sédiments et de sismique Haute Résolution[END_REF].

According to the Hamilton measurements, the variation of absorption with frequency is a linear relationship. In contrast, according to Kibblewhite a linear relationship can not be established over the entire frequency range found in a seismic application to ultrasound [START_REF] Kibblewhite | Attenuation of sound in marine sediments : a review with emphasis on a new low-frequency data[END_REF]. This would be due, in part, to the dissipation caused by the fluid viscosity. These conclusions are consistent with the Biot theory but there is still a debate about the frequency dependence of absorption coefficient in the underwater acoustics community.

Characterization of seabed by acoustic methods

Why are the acoustic methods used ?

The knowledge of the seabed is an essential problem in ocean exploration. For a long time, the investigations were based on direct observations in the medium or by sampling. The use of underwater acoustics has improved these technologies because the acoustic wave has a high capacity in propagating in water at long distances. Thus, it makes it easy to convey information, especially the information useful for measuring the characteristics of the various constituents of the marine environment, particularly the shape and the nature of the bottom [START_REF] Lurton | Présentation et applications[END_REF]. This is why many technological and acoustic means have been developed to characterize the seabed. These acoustic methods are quite recent since the first reliable practical achievements dates back to the 20th century. The work presented in this thesis is a part of those methods and aims to develop an approach for geoacoustic inversion.

CHAPTER 1. SEDIMENTS AND ACOUSTICS

In this part of the chapter, we propose to begin by introducing the main sonar systems used to identify the sedimentary structures as well as their instrumentation and their algorithms. Then we present the methods of geoacoustic inversion.

Geoacoustic systems

The acoustic measurement systems for seabed characteristics are predominantly of active type (controlled emission). There are several acoustic devices which have been developed in the field, such as the echo sounder, the multibeam echosounder, the sub bottom profiler, the side scan sonar, etc. We can classify these systems according to the intended objective (frequency) [START_REF] Rakotonarivo | Caractérisation quantitative des fonds marins à l'aide des sondeurs de sédiment (1 kHz -10 kHz) : modélisation directe pour l'inversion sur données réelles[END_REF] or according to the method used (mode) [START_REF] Viala | Inversion géoacoustique temps réel de signaux large bande par grands fonds[END_REF].

Most geaoacoustic measurements are made in remote sensing which refers to the extent or acquisition of information about an object through a measurement instrument having no contact with the object studied. This is made possible by the fact that the studied objects (seabed or sediments) reflect and scatter waves at different wavelengths and intensities according to their nature. The analysis of the signal echo sent after its interaction with the bottom makes it possible to extract information about the geological and the geoacoustic characteristics of the seabed, such as the depth, the density, the porosity, the roughness or the sound speed.

A first step of this algorithm consists in processing the measured temporal echo to extract the observed data. The characteristics of the seabed are parameters that modify the observed data. Finding these parameters allows to find the inverse transfer function associate with the measurement (Figure 1.9).

Frequency influence

The acoustic propagation features depend on the frequency. Information about the hardness or the roughness of the bottom obtained by a 200 kHz measurement is not necessarily valid at 2 kHz. Knowing the characteristics of the seabed is thus related to the frequency of use. For seabed exploration, two major classes of devices can be distinguished :

-Medium and high frequency sonar systems (>10 kHz) which map the surface of the seabed and have a fair estimate of the undersea feature. -The low frequency sonar (<10 kHz) which provides the seabed profiles. The source used is generally impulsive with high power and low frequency (LF) so that the emitted waves can penetrate into the buried layers of the seabed just a few meters for the sediment sounder and a few kilometres deep for the seismic.

CHARACTERIZATION OF SEABED BY ACOUSTIC METHODS

21

In seismic the instrumentation includes an explosive source (Boomer) or an electrical discharge (high resolution sparker) in emission that sends low frequency signals (100 Hz or below). The receiving antenna is an hydrophone array towed at the back of a boat, varying in length from 100 m to 4 km. The measurement configuration is essentially bistatic (see in the next paragraph). When the incident angles are less than the critical angle, the signal measured by the reception is the subject of seismic reflection. Beyond the critical value, the acoustic field is fully transmitted in the structure. The wave travels parallel to the surface, then it is sent back to the receiving sensor. This process is the subject of seismic refraction (Figure 1.10). . The sediment echo sounder consists of electro acoustic transducers which work alternatively in emission and then reception. The signals issued generally have frequencies in the range 1 kHz-10 kHz, and penetrate to the vertical relative to the ground and are at the same position. The received signal corresponds to backscatter, i.e the energy returns in the emission direction.

Geometry influence

In remote sensing, the sonar systems emit signals which are reflected and scattered by the bottom, sent back to the surface and recorded by one or numerous hydrophones. If the system is placed near the sound source or in the same place, we speak of a monostatic system (Figure 1.11) such as the echo sounder, the multibeam sounder and the side scan sonar, otherwise it is a bistatic system such as seismic (Figure 1.12).

Monostatic configuration

The monostatic system has an interest in applications such as dredging, finding the sediment thickness or characterizing the seabed (stratification, slope, features). An example of a monostatic system is the echo sounder. It has certain advantages, namely, the characteristics of the emission antennas are well known, the directivity can be controlled, the offset between the sensors is zero. The frequency range of the sediment sounders (usually 1 kHz-10 kHz) is higher than in seismic which means that the resolution is better but the signal penetrates less into the background : few kilometers in conventional seismic and a few tens of meters for sediment sounders. Moreover, from a single shot, they can be traced back only to 1-D information while seismic devices provide access to 2-D information.

One of the disadvantage of the monastic system, in the context of seafloor characterization, is the vertical incident. It might cause confusion in estimating geoacoustic parameters, e.g between sound speed and thickness of layers in the stratified medium (Echo sounder) or give access only to the scattering signals (Multibeam sounder, Side scan sonar) which is often a too poor data to perform efficient geoacoustic inversion. In fact, the information gathered by monostatic methods are very general and unsuitable for specific use such as the transmission sonar system. So, if we want access to the parameters of the stratified medium, reflected signals measurement must be done with bistatic configuration.

Bistatic configuration

The characteristics of bistatic configuration is the use of one or more spatially separated receivers. This ensures the exchange of track data and the appropriate fusion of all the available data in order to achieve the optimal detection performance (Figure 1.12). The bistatic configuration provides some well known advantages compared with monostatic configuration like enlarged detection areas and increased probability of detection. More precisely, in geoacoustic inversion, the advantage is the use of a different detection angle and different propagation waves in a stratified medium.

Geoacoustic inversion methods (a) Definition

The inverse problem consists in estimating the parameters describing a physical system, given the observed values of this system at certain measurement points. Direct access to the various parameters describing the physical system can be difficult, sometimes impossible. Nevertheless there are often physical quantities related to those parameters that are more easily accessible by measurement. In underwater acoustics, these physical quantities might be the sound pressure, the backscattered field, the reflection coefficient, etc. The parameters could be those describing the water column (bathycelerimetric profiles, water depth), the seafloor, or the localization of a target.

To model this type of problem, we assume that the parameters to be estimated are the components of a vector m. The components of d obs = [d 1 , d 2 , . . . , d N ] are the data that are observed. The relationship between the observations d obs and the parameters m is assumed to be functional :

d obs = G(m).
(1.10)

This relation implies that the functional G is perfect. But, in real world, there is some discrepancies between the data and the predictions of G that can come from the measurements (noise, problems with the sensors,. . . ) or from the model itself (approximations, hypotheses,. . . ). These discrepancies might be modeled by adding noise to the previous functional :

d obs = G(m) + n. (1.11)
If the inverse problem is linear, the vector parameter m is simply estimated by finding the inverse operator G -1 :

m = G -1 (d obs ).
(1.12) But unfortunately, this is rarely the case and most of the inverse problems in underwater acoustics are non-linear.

From a general point of view, an inverse method is composed by the following steps :

1. choice of the parameters to be estimated and a priori on these parameters, CHAPTER 1. SEDIMENTS AND ACOUSTICS 2. selection of the measurements and the data to be inverted with relations to the parameters, 3. definition and computation of the forward model G, 4. choice of a cost function to measure the difference between the data and the predictions, 5. choice of an optimization strategy to minimize the cost function.

In the following, these points are briefly described.

(b) Choice of the parameters

Obviously, the inverse problem is firstly conditioned by the parameters we would like to estimate. The whole strategy is different if, for instance, we want to localize a whale from its acoustic emission or if we are interested by the micro-roughness of a seafloor. The definition of parameters affects the main characteristics of the measurements : frequency, geometry acquisition (monostatic vs bistatic), signal shape, passive or active, etc.

In this initial step of the inverse problem, it is possible to introduce a priori on the parameters. This comes from the knowledge of the problem physics. Two different types of a priori might be introduced :

-bounds on the parameters values, -relationships between different parameters.

(c) Measurements and data

The measurements made to produce the data d obs to be inverted must be closely linked with the choice of the parameters. But, in reality, two different cases may exist. First, a specific campaign can be designed and realized specially to produce the data to invert with close relation to the parameters. Second, one can used data already acquired by a campaign and test their match with a new problem.

In order to test the developed algorithms, simulated data are very often use. They offer the advantage of producing data from a perfectly known environment and therefore to check the performance of the inverse method.

(d) Forward model

The objective of the forward (or direct) model is to produce simulated data d sim which must be sensitive to the parameters and which are compared to the observed data. Depending on the characteristics of the experiments, different model can be used to this purpose. In underwater acoustics, a lot of numerical codes are available to solve the wave equation in order to predict the pressure field in the water column. These codes are very well described and documented in the reference textbook [START_REF] Jensen | Computational ocean acoustics[END_REF].

Essentially, these codes can be grouped in five families : fast field programs (FFP), normal modes (NM), rays, parabolic equation (PE), and direct finite-difference (FD) or finite-elements (FE). The three last families are naturally adapted to range-dependent environments, i.e. to environments which characteristics vary with range. However, it is possible to extent FFP or NM codes with coupled solutions to treat range dependent problems.

Among all the characteristics of the experiments to be taken into account for the choice of the forward model, the frequency is certainly the most important. Regular choices are often discriminate by :

-in the case of few kilohertz or above rays and FFP are the most used, -in the case of a kilohertz or below NM, PE, and FD/FE are predominant. These codes are now well known, with their performances and limits. But they might not be adapted to every measurements. For instance, a multibeam echosounder records the acoustic signal backscattered by a rough and heterogeneous seafloor. If we want to build an inversion method on this type of data, a specific forward model must be used.

(e) Cost function

Once the experiments have been made, producing d obs , and the forward model chosen and run, producing d sim with parameters m, the next step is to compare the observation and the predictions in order to obtain the best parameters vector possible. This operation is usually performed through the definition of a cost function f (e) where e is the prediction error : e = d obsd sim .

(1.13)

This function might take different forms, usually grouped in two families.

The first family is based on the definition of the norm of the prediction error. In discrete case, this norm can be expressed in different ways :

-L p norm : f (e) = e p = [ e i p ] 1/p , -weighted L p norm : f (e) = e w,p = [ w i × e i p ] 1/p where w is a weighting window, -L ∞ norm : f (e) = e ∞ = max i e i . The choice of the norm depends on the problem and on the data : the higher order the norm is, the more erroneous data are taken into account.

The second family is based on the correlation between predictions and observations. Different functions can be defined, here again depending on the problem and the data. Among, all the possible definitions, the most common are :

-simple correlation : f (d sim , d obs ) = d obs * • d sim d obs 2 • d sim 2 , -Bartlett correlator : f (d sim , d obs ) = d sim * • C • d sim with C = d obs • d obs * , -maximum likelihood : f (d sim , d obs ) = 1 (d sim ) * • C -1 • d sim .
In these expressions, the exponent * indicates the complex conjugaison.

(f ) Optimization

The solution to the inverse problem, or the best parameters vector m, correspond to an extremum of the cost function : a minimum for the ones based on norm, and a maximum for the ones based on correlation. Therefore, an optimization strategy is needed. This final step is crucial and might be difficult to resolve because many inverse problems are ill-conditioned. The most common issue is when the cost function have local extrema. Therefore, the optimization strategy must be able to get rid of theses local extrema to converge to the global extremum. Among all the optimization strategies used in the underwater acoustics community, the following are the most common :

-gradients (simple or conjugate), -Newton, -genetic algorithm [START_REF] Gerstoft | Ocean acoustic inversion with estimation of a posteriori probability distributions[END_REF], -simulated annealing [START_REF] Kuperman | Non linear inversion for ocean bottom properties[END_REF]. The choice of one optimization algorithm depends on the nature of the problem : quasilinear or strongly non-linear, number of parameters, structure of the cost function,. . . Another approach to this problem is the probabilistic inversion, known as the Bayesian methods. In this case, the objective is not to obtain the best estimation of the vector parameter m, but to estimate its posterior probability density (PPD). This approach gives in fine an idea of the "quality" of the parameters estimation and takes into account the data and model uncertainties [START_REF] Tarantola | Inverse Problem Theory And Methods For Model Parameter Estimation[END_REF].

The objective is to obtain the PPD expressed in Bayes' theorem :

P (m|d obs ) = P (d obs |m)P (m) P (d obs ) . (1.14)
Once the experiment has been carried out and the data obtained, P (d obs ) is a constant factor which can be absorbed into the normalization, and the conditional probability density function (PDF) P (d obs |m) can be interpreted as a function of m known as the likelihood function. With some assumptions on the data and the model, the likelihood function has an analytic expression and can be sampled, for example with Metropolis-Hastings algorithm to provide the desired PPD.

The image source method

One major drawback of the inversion methods described in the previous section is that they are often very intense in computation. Compared to these approaches, the inverse geoacoustic method studied in this thesis, the image source method (ISM), initially developed by S. Pinson, has a very low computational cost [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF]. Before describing it in details in the following chapters, it is now replaced within the general context of inverse problems described in the previous section.

The parameters that are sought by the original ISM are : the number of layers of the superficial seafloor, the thicknesses of these layers, and the sound speed of P waves inside these layers. In this thesis, we propose to add two other parameters estimation to ISM : the density of the layers and the S waves sound speed. This is described in chapter 6. This choice of parameters come from the description of the seafloor made in sections 1.3 and 1.4 and is the most common choice for the geoacoustic seabed inversion.

To estimate these parameters, the data to invert must be sensitive to them. Based 1.6. CONCLUSION on the remarks made in section 1.5.2 about monostatic and bistatic systems, it appears that the bistatic configuration is the most able to provide such data. Indeed, to discriminate thicknesses from sound speed effects on travel times, it is mandatory to have some refraction on the acquisition system. Therefore, the ISM is based on the analysis of acoustic signals reflected by the seafloor in the time domain.

The frequency used must be low enough to investigate the first meters of sediments but high enough to have a satisfactory resolution. With respect to this inherent compromise between range and resolution, we decided to work in the mid-frequency range, i.e. between 1 and 10 kHz. This frequency band corresponds to wavelength between 15 cm and 1.5 m and allows a sufficient signal penetration in the seafloor.

The main difference between ISM and other geoacoustic inversion methods arises from its forward modeling. Other geoacoustic inversion methods based on the same data uses either all the travel times between the source and the hydrophones array or computes the reflection coefficient in frequency domain. The forward modeling of ISM implicates a drastic reduction in the number of data to invert and therefore, to a very light inversion scheme. In fact, no obvious cost function and optimization strategy needs to be defined. This is explained in details in section 3, 4, and 5 of this thesis.

Conclusion

To analyze the marine sediments by acoustic methods, it is necessary to have a good understanding of the direct problem, i.e how the waves propagate in the known environment. To this end, in this first chapter we have described the marine sediments that represent the environment to be studied, the acoustic means used and its inversion algorithms.

Chapter 2

Synthetic and real signals

Introduction

As described in section 1.5.3, the definition of an efficient geoacoustic inversion scheme is strongly dependent on the measurements made and on the type of data provided. The ISM studied in this thesis is based on the analysis of acoustic signals reflected by a layered seafloor in the mid-frequency band and analyzed in the time domain. Before describing the ISM, this chapter is devoted to the presentation of the data used in this work. In fact, two types of data are presented.

First, synthetics data are computed. The objective here is to produce data that mimic real ones in order to check the efficiency and the quality of ISM. Different types of acquisition geometries and seafloors are used to produce these data.

Second, two real dataset are used : one acquired in a water tank with artifical seafloor, and one coming from at-sea data acquired by NURC (Nato Undersea Research Center, La Spezzia, Italy) in 1998.

CHAPTER 2. SYNTHETIC AND REAL SIGNALS

Synthetic signals

The signals to synthesize are the reflection of a broadband signal by a stratified seafloor medium. There are several approaches (see section 2.2.1) to calculate the impulse response of such medium. The signals generated by the direct model should be the closest possible to real signals and sensitive to the parameters which we search to estimate. For this, the choice of the direct model, the source type, the model of the bottom have to be taken into account accurately.

In practice, the measurements are made with a source whose dimensions are limited. If these dimensions are small compared to the wavelength, one can make the assumption of a point source emitting a spherical wave. The problem with this approximation is that the plane interface on which the wave is reflected is geometrically different from spherical waves. The decomposition of plane waves, with Sommerfeld integral, into an infinite and continuous sum of plane waves solves the problem. Then, we found the reflection coefficient of plan wave in stratified medium with different approaches presented in section 2.2.2 for each case of layers type, fluid or solid, under the Born approximation.

To explain this approach in details, we first recall the propagation equations of spherical waves in a free field. Then, we present the pressure field recorded after a reflection on a fluid interface in a first time and then at a solid interface, recalling each time the values of the reflection coefficient to which it corresponds in a stratified medium.

Reflection of an acoustic spherical wave on a layered seafloor

Acoustic wave in free field

The combination of the linearized equations (mass conservation and Euler) and the state equation leads to the wave equation for pressure p :

∆p - 1 c 2 ∂ 2 p ∂t 2 = 0 , (2.1)
where c is the sound speed in the medium. In the spherical coordinates system (r, θ and ϕ), the Laplacian is written as :

∆F = 2 r ∂F ∂r + ∂ 2 F ∂r 2 + 1 r 2 sin ϕ ∂ ∂φ sin ϕ ∂F ∂φ + 1 r 2 sin 2 ϕ ∂ 2 F ∂θ 2 .
(2.2)

In the case of a spherical wave, the radiation is isotropic and the laplacien is reduced to :

∆F = 2 r ∂F ∂r + ∂ 2 F ∂r 2 = 1 r ∂ 2 (rF ) ∂r 2 . (2.
3)

The wave equation can therefore be written as :

∂ 2 (rp) ∂r 2 - 1 c 2 ∂ 2 rp ∂t 2 = 0 . (2.4)
So the solution for a progressive spherical harmonic wave towards the increasing r is written as :

p = A r e i(wt-kr) , (2.5)
where ω is the pulsation and k is the wavenumber. Now, let us consider a spherical source at the origin of cartesian system of coordinates, the pressure in the plane z = 0, omitting the term e -iωt and an amplitude factor is as follows :

p(r) = e ik r r , (2.6)
where r is the distance x 2 + y 2 .

According to [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF] the extension of the pressure field in the double Fourier integral is described in the (x, y) plane as follows :

e ikr r = +∞ -∞ A(k x k y )e (i(kxx+kyy)) dk x dk y , (2.7)
Finally, the spherical wave can be decomposed into an infinite sum of plane waves [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF], the pressure field in the case of a positive and negative z is written as follows :

For z ≥ 0 exp ikD D = i(2π) -1 +∞ -∞ k -1 z exp(i(k x x + k y y + k z z))dk x dk y , (2.8)
For z ≤ 0

exp ikD D = i(2π) -1 +∞ -∞ k -1 z exp(i(k x x + k y y -k z z))dk x dk y , (2.9) with D = x 2 + y 2 + z 2 .

Acoustic wave in the case of a plane interface

The development of the expression of a spherical wave in the free field on a plane wave basis addresses the problem of the different symmetry between this spherical wave and the plane interface. Indeed, the influence of the medium below the interface on the reflected field is then simply taken into account by multiplying the previous integrand by the reflection coefficient R(k z ) of plane wave, where k z is the vertical projection of the wave vector. The details of calculation of this equation are given in reference [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF].

We assume that the source at (x = 0, y = 0, z = z 0 ) and the plan z = h s is the first interface of the layered medium. Using the decomposition of the spherical wave into continuous infinite sum of plane waves in space, the phases of the latter in the integral has become after reflection on the interface k x x + k y y -k z (z + z 0 ). Hence the recorded pressure field is given by :

p r = i(2π) -1 +∞ -∞ k -1 z R(k z )exp(i(k x x + k y y + k z (z + z 0 )))dk x dk y (2.10) CHAPTER 2.

SYNTHETIC AND REAL SIGNALS

This equation can be simplified by using polar coordinates :

k x = ξ cos ψ, k y = ξ sin ψ, ξ = k 2 x + k 2 y ,
x = r cos ϕ, z = r sin ϕ, dx dy = rdrdϕ .

(2.11)

The equation 2.10 becomes :

p r = i(2π) -1 +∞ -∞ k -1 z R(k z )exp(ik z (z + z 0 ))ξdξ 2π 0 exp(iξr cos(ϕ -ψ))dψ. (2.
12)

The part of integral in ψ can be developed using the Bessel and the Hankel functions. After calculation, the reflected pressure is given by :

p r = i(2) -1 +∞ -∞ k -1 z R(k z )H 1 0 (ξr) exp(ik z (z + z 0 ))ξdξ, (2.13)
where H 1 0 is the Hankel function of the first kind. This formula is known as the Sommerfeld integral. After changing the integration variable, it may also be written as :

p r = ik π 2 -i∞ 0 J 0 (kr sin(θ))R(θ)e (ikz(z+z 0 ) cos(θ)) sin(θ)dθ, (2.14)
where J 0 is the Bessel function of the first kind. This formulation has two important advantages : -First, R(k z ) is the reflection coefficient of a plane wave which can be easily computed for different configurations of seabed even a complex one (section 2.2.2) and with different approaches. -Second, the integral boundary contain an imaginary part which corresponds to inhomogeneous waves. Therefore, complex phenomena such as head waves are modeled with this approach.

This formula is widely used in this thesis because all the simulated signals are computed trough a numerical evaluation of the Sommerfeld integral by using equation 2.14 to find the transfer functions [START_REF] Holland | Resolving mesoscale seabed variability using reflection measurements from an autonomous underwater vehicle[END_REF].

Born approximation of Sommerfeld integral in a fluid stratified medium

The Sommerfeld integral (equation 2.14) describes the reflection of a spherical wave by a fluid stratified medium. In the case of small impedance difference between the layers and of incident angle under the critical one, the Born approximation can be used [START_REF] Clayton | A born-wkbj inversion method for acoustic reflection data[END_REF][START_REF] Rajan | A comparison between the born and rytov approximations for inverse backscattering problem[END_REF]. Sommerfeld integral can be simplified by the fact that we no longer take into account the multiple reflection. This assumption is valid because the value of amplitude in the reflection is usually higher by about 2 magnitude orders in the first term of multiple reflections [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF]. To check the validity of this approximation, the analytical expression of the combination of the Sommerfeld integral and the Born approximation in a stratified medium are developed in details in reference [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF]. In this work, under the above mentioned hypotheses, the Sommerfeld integral is evaluated by the method of stationary phase and the expression can be rewritten according to θ

(n) q = sin -1 c q c 0 sin(θ (n) 0 ) which
is the incident angle in the layer q and corresponding to the reflection of the wave on the interface n as follows :

p(r 0 , r, ω) = N -1 n=0 R θ n 0 1 x cos 2 θ n 0 c 0 sin n 0 n q=0 2c q h q (cos 3 θ n q ) × (2.15) exp[i(k 0 sin (n) 0 x + n q=0 2k 0 c 0 c q cos θ n q h q )] (2.16)
where R is the amplitude factor which includes the local reflection coefficient, h is the thickness layer, c is the sound speed, k is the wavenumber and N is the number of layers. With this formulation, the only element that can not be obtained directly is the value of incident angles. This angle can be obtained by Fermat principle :

θ (n) 0 = argmin θ (n) 0 (k 0 sin (n) 0 x + n q=0 2k 0 c 0 c q cos θ (n) q h q ).
(2.17)

Consequently, in this configuration of relatively low impedance contrasts between the layers and with incident angles far from the critical angle, the multiple reflections inside the sediment stack may be neglected. Figure 2.1 shows the rays considered, in the case of Born approximation, in the example of a 2 layers over a basement system in which only 3 direct reflections on the three interfaces are included.

Water (c 0 ) Layer 1 (c 1 , ρ 1 ) Layer 2 (c 2 , ρ 2 ) Basement x z Source (x s = 0, z s = 0) Receiver z K x K h s θ (0) 0 θ (0) 0 h r = h s -z K h s + h 1 θ (1) 1 θ (1) 0 θ (1) 0 h 1 h s + h 1 + h 2 θ (2) 2 θ (2) 0 θ (2) 0 h 2 Figure 2
.1 -Reflection of a wave on a stratified medium in which only the first reflections on the interfaces are shown.

Plane wave reflection coefficient of a complex seafloor

The case of a simple interface

Let us consider two fluid semi-infinite media with densities ρ 0 and ρ 1 and sound speeds c 0 and c 1 separated by a plane and horizontal interface at z = h s (Figure 1.7). The reflection and transmission coefficients (equation 1.8) can be written as :

R = ρ 1 c 1 cos(θ (0) 1 ) - ρ 0 c 0 cos(θ (0) 0 ) ρ 1 c 1 cos(θ (0) 1 ) + ρ 0 c 0 cos(θ (0) 0 ) , (2.18) T = 2ρ 1 c 1 cos(θ (0) 1 ) ρ 1 c 1 cos(θ (0) 1 ) + ρ 0 c 0 cos(θ (0) 0 ) , (2.19) θ (0)
0 is the incident angle at the medium 1/medium 2 interface and θ

(0)
1 is the transmission angle in the medium 2.

The case of a simple interface in fluid-solid

If the medium below the interface is solid, it supports the propagation of two types of waves : shear wave (sound speed c s ) and longitudinal wave (sound speed c p ). In this case, the reflection coefficient can be written [START_REF] Viala | Inversion géoacoustique temps réel de signaux large bande par grands fonds[END_REF] :

R = (Z T -Z 1 ) (Z T + Z 1 ) , (2.20)
where Z T is the impedance of the second media :

Z T = Z p sin(2θ p ) 2 + Z s sin(2θ s ) 2 .
(2.21) θ p and θ s check Descartes law :

k 1 sin(θ (0) 0 ) = k s sin(θ s ) = k p sin(θ p ) . (2.22)

Reflection coefficient in the case of stratified fluid medium

To simplify the problem, we start this paragraph by the computation of the reflection coefficient of one layer above a semi-infinite basement (Figure 2.2). Different methods exist to find it. One of these methods is presented in this section [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF]. The successive reflection and transmission coefficients at each passage by interface and their phase shift are used to this purpose. The reflected field is written as the sum of reflections and multiple reflections in the layer :

Water (c 0 , ρ 0 ) Layer 1 (c 1 , ρ 1 ) Basement (c 2 , ρ 2 ) x z Source (x s = 0, z s = 0) T 01 T 01 R 12 T 01 R 2 12 R 10 R 01 T 01 R 12 T 10 T 01 R 2 12 R 10 T 10 T 01 T 12 T 01 R 12 R 10 T 12 h 1 h 2 θ 1 θ 1 θ 1 θ 2 θ 1 θ 2 θ 3
R = R 01 + T 01 T 10 R 12 exp(2iφ 1 ) + T 01 T 10 R 10 R 2 12 exp(4iφ 1 ) + T 01 T 10 R 2 10 R 3 12 exp(6iφ 1 ) + ... = R 01 + T 01 T 10 R 12 exp(2iφ 1 ) × ∞ n=0 (R 12 R 21 exp(2iφ 1 )) n ,
(2.23) where R nm is the reflection coefficient from medium n to medium m and T nm is the corresponding transmission coefficient, φ 1 = k 1 h 1 cos(θ 1 ), k 1 is the wavenumber in the first layer and h 1 is the thickness of the layer. By using geometric series, the relation 1 + R = T and the relation R mn = -R nm , the reflected and transmission coefficients field are written as :

R = R 01 + R 12 e (2iφ 1 ) 1 + R 12 R 01 e (2iφ 1 ) = Z 1 (Z 2 -Z 0 ) -i(Z 2 1 -Z 0 Z 2 ) tan(φ 1 ) Z 1 (Z 2 + Z 0 ) -i(Z 2 1 -Z 0 Z 2 ) tan(φ 1 ) , (2.24) T = T 01 -T 12 e (2iφ 1 ) 1 + R 12 R 01 e (2iφ 1 ) .
(2.25)

Reflection coefficient in the case of multiple fluid layers

Let generalize the previous results to the case of medium with n -1 fluid layer over a semi-infinite basement. For this case, two different approaches to find the total reflection coefficient are presented.

The first method [START_REF] Jensen | Computational ocean acoustics[END_REF] is based on the equation 2.24 in which we assume that the calculation of the deepest reflection coefficient of the three last layers is :

R (n-2)n = R (n-2)(n-1) + R (n-1)n e (2iφ n-1 ) 1 + R (n-2)(n-1) R (n-1)n e (2iφ n-1 ) . (2.26) CHAPTER 2.

SYNTHETIC AND REAL SIGNALS

Knowing R (n-2)n , we compute the next reflection coefficient as follows :

R (n-3)n = R (n-3)(n-2) + R (n-2)n e (2iφ n-2 ) 1 + R (n-3)(n-2) R (n-2)n e (2iφ n-2 ) . (2.27)
This process is iterative and finally the total reflection coefficient field is obtained by :

R 0n = R 01 + R 1n e (2iφ 2 ) 1 + R 01 R 1n e (2iφ 2 ) .
(2.28)

To understand better this algorithm, we take the example of a stratified medium with n = 4 (figure 2.3).

Using equation 2.24 we find the total reflection coefficient of the three last layers R 24 as follows :

R 24 = R 23 + R 34 e (2iφ 3 ) 1 + R 23 R 34 e (2iφ 3 ) .
(2.29)

Having found R 24 , the reflection coefficient R 14 is given by the simple formula :

R 14 = R 12 + R 24 e (2iφ 2 ) 1 + R 12 R 24 e (2iφ 2 ) , (2.30)
The total reflection coefficient of the whole stratified medium is found as follows :

R 04 = R 01 + R 14 e (2iφ 1 ) 1 + R 01 R 14 e (2iφ 1 ) , (2.31)
where R 01 is a reflection coefficient at the water/layer 1 interface and R 14 is reflection coefficient of all the layers below.

The second method to calculate the total reflection coefficient of a plane wave over a stratified medium with an arbitrary number of fluid layer n is based on impedances [START_REF] Brekhovskikh | Fundamentals of Ocean Acoustics[END_REF]. To explain briefly the principle, let consider that a plane wave is incident on the water/layer 1 interface with an angle θ 0 (Figure 2.4). The main object is to determine the global reflection coefficient. So, thanks to the continuity of normal impedances in each layer it is possible to "bring the interface (n-1/n) back to interface 0" by using recurrence formulas to obtain finally the equivalent input impedance Z

(1) in of this multilayered system. For any layer p, with p ∈ [1, n -1], we have :

Z (p) in = Z p Z (p+1) in -iZ p tan(φ p ) Z p -iZ (p+1) in tan(φ p ) , (2.32)
The recurrence is initialized with

φ n = k zn h n , h n = 0 and Z (n) in = Z n .
Consequently, the incidental wave sees only a semi-infinite medium characterised by impedance Z

(1) in and the reflection coefficient is simply given by : R = Z

(1)

in -Z 0 Z (1) in + Z 0 .
(2.33)

Water (c 0 , ρ 0 ) Layer 1 (c 1 , ρ 1 ) Layer 2 (c 2 , ρ 2 ) Layer 3 (c 3 , ρ 3 ) Basement (c 4 , ρ 4 ) R 01 R 12 R 23 R 34 R 24 R 14 R 04 h 1 h 2 h 3 h 4 Figure 2.3 -Reflection coefficient in the case of multiple fluid interfaces water (c 0 , ρ 0 ) Layer 1 (c 1 , ρ 1 ) Layer 2 (c 2 , ρ 2 ) Layer 3 (c 3 , ρ 3 ) Layer n -2 (c n-2 , ρ n-2 ) Layer n -1 (c n-1 , ρ n-1 ) Basement n (c n , ρ n ) Z (0) in Z (1) in Z (2) in . . . Z (n-2) in Z (n-1) in θ 0 θ n h 0 h 1 h 2 h n-2 h n-1
Figure 2.4 -Parameters for the calculation of the reflection coefficient of a system of layers.

Other aspects

Presence of elastic layers

When we simulate the presence of solid layers within the sediment stack, we carry out the numerical evaluation of equation 2.14 with global reflection coefficient R calculated for this kind of bottom. To perform this, we could use a specific software such as OASES. However, we prefer to utilise an algorithm based on a transfer matrix [START_REF] Cervenka | A new efficient algorithm to compute the exact reflection and transmission factors for plane waves in layered absorbing media (liquids and solids)[END_REF]. The fact that acoustic waves propagate in each layer is revealed through a matrix that relates the velocity potential for both compressional and shear waves at the input with their values at the output. The limit conditions at each interface are also expressed in a matrix formalism that ultimately brings about a compact formulation appropriate for a reliable calculation.

Presence of dissipation layers

The absorption is taken into account by two different ways : making a complex sound speed or a complex wave number. According to the general chosen agreement in underwater acoustic, a harmonic plane wave is written as follows :

p = P 0 e i(kr-ωt) , (2.34)
where P 0 is the amplitude, k is the wave number and ω is the pulsation. The absorption results in the fact that the wave number is then complex :

k = k r + iN e , (2.35)
where N e is the attenuation in Np/m.

The absorption can also be described by taking the complex sound speed :

c = c r + ic i , (2.36)
with c i > 0. Since c i c r , we can write :

k ω c r + i ωc i c 2 r where c i = αc 2 r ln(10) 40πf , (2.37)
where f is the frequency and α the attenuation in dB/m. In some publications, the "loss parameter" is defined as :

δ = c i c r .
(2.38)

Geometric configuration used

The ISM has already been tested on real at-sea data with two different geometric configurations : a vertical antenna (see references [START_REF] Pinson | Sound speed profile characterization by the image source method[END_REF][START_REF] Drira | Image source detection for geoacoustic inversion by the teager-kaiser energy operator[END_REF]) and a horizontal antenna (see references [START_REF] Pinson | Range dependent sediment sound speed profile measurements using the image source method[END_REF]). The aim of this part is to reproduce these two configurations with synthetic signals in order to understand the effects of each parameter variation on the ISM. So, we can analyse the parameter variation according to the emitted signals, the wave distance propagation, the type of antenna and the effect of different type of geoacoustic configurations.

Configuration 1 : vertical array of hydrophones.

The simulation is performed in order to present the Scarab 98 [START_REF] Holland | High resolution geoacoustic inversion in shallow water : A joint time and frequency domain technique[END_REF] configuration, a real experiment presented in the section 2.3.2, with a vertical antenna of 15 hydrophones spaced irregularly and a length of 64 m. The hydrophone coordinates are presented in the table 2.1. The source is located at 150 m over the bottom and 200 m from the array of hydrophones (Figure 2.5).

Configuration 2 : horizontal array of hydrophones

The simulation is performed in order to present the Clutter 09 configuration. In spite of the fact that the real experiment is not presented in this thesis, the choice of this configuration is used to test the efficiency of the ISM algorithm to estimate the density, the sound speed c s of shear wave and the absorption coefficients in a simple medium in chapter 6, and to test the effect of the horizontal antenna on the new ISM algorithm. This configuration has an other advantage : thanks to this device, an inversion of range-dependent medium can be done.

The geometric configuration is as follows : a broadband source at 12 m over the CHAPTER 2. SYNTHETIC AND REAL SIGNALS seafloor, the array is horizontal made of 15 hydrophones at the source height as the source and 24 m away from it. The hydrophones are spaced regularly, 1 m between each two adjacent hydrophones so its length is 14 m from (x = 24, z = 0) to (x = 38, z = 0) (Figure 2

Geoacoustic configuration used

In this section, different types of seafloor are used to simulate bottoms with parameters variation : for each type of seafloor we give the sound speed, the density, the thickness and the absorption coefficient. The evolution of sound speed according to depth is then called SSP (Figure 2.7). The aim of the different type of the seafloor is to understand the impact of these parameter changes on ISM inversion results in the next chapters.

Seabed 1

The seafloor is composed of 9 fluid layers on a basement, the geoacoustic parameters of these layers are presented in the table 2.2. The absorption is assumed zero in the medium. The choice of a large number of layers is to test the efficiency of the algorithm of the ISM with depth, the effect of the number of layers in the resolution of inversion and finally to study the effect and the propagation of an error detection on the first layers to the rest of inversion parameters (Chapter 3).

An example of the SSP of the seabed 1 is presented in the figure 2

Seabed 2

The seafloor is composed of two fluid layers over a basement in the first case (then simply called seabed 2) and two solid layers over a basement in the second case (called seabed 2 solid). The geoacoustic parameters are presented in table 2.3. The absorption coefficient in the layers are taken into account in this example. The choice of a small number of layers is to test the accuracy of the algorithm of the ISM with a simple seabed, to simplify the algorithm of inversion in chapter 6 to find the shear wave, the density and the absorption coefficients of layers.

CHAPTER 2. SYNTHETIC AND REAL SIGNALS

Seabed 4

The seafloor is composed of 4 fluid layers on a basement, the geoacoustic parameters of these layers are presented in the table 2.5. The absorption is assumed zero in the medium.

The choice of the geoacoustic configuration 3 and 4 is to test the effect the error model of the algorithm of the ISM with different seabed in Chapter 3.

Layers Sound speed (m/s) Thickness (m) Water 1500 . layer 1 1540 6.5 layer 2 1600 4 layer 3 1700 8 layer 4 1720 7

Table 2.5 -The geoacoustic parameters of the layered medium labeled seabed 4. The absorption coefficient of layers are assumed zero.

Examples of signals simulated

The signal recorded in the hydrophone K is composed of the reflected wave on the stratified medium and the direct wave between the source and the receiver is given by.

s K (t) = G 0 (t) * f (t) + H K (t) * f (t),
(2.39)

where f (t) is the emitted signal and G 0 (t) is the Green function in a homogeneous water sound speed c 0 . H K (t) is the transfer functions of acoustic wave reflections on the layered medium which is found by using the numerical evaluation of the Sommerfeld integral equation 2.14.

In the following, different examples of such simulations are presented with different configurations and seabeds.

The signal emitted (Figure 2.8) in the case of configuration 1 and seabed 1 is the real signal recorded in the first hydrophone. Equation 2.39 is used to compute the synthetic signals with the configuration 1 and seabed 1. Though the absorption is assumed zero in the medium, we find that the amplitudes of the multiple reflections are too weak in the simulated signal (Figure 2.9) to be detected compared to the first reflections. This is probably due on the one hand to the low impedance contrasts between the interface and on the other hand we have almost a total reflection. This result confirms the remark made on the Born approximation of the Sommerfeld integral subsection 2.2.1.

Real signals

Tank experiment

Tank description

To have real signals in a controlled environment with a precise knowledge of the ground truth, we have realized an experiment in a tank at ISEN (Institut Supérieur de l'Electronique et du Numérique, Brest) [DMB + 15].

The experimental tank used in this study is shown in figure 2.12. As depicted in figures 2.13 and 2.14, it consists of a basin that is 3 meters long and 2 meters wide with a water depth of 1 meter. The tank bottom is made of canvas rendered by soft plastic. The tank is equipped with 2 metal arms immersed in water and carrying sensors with an handle. The displacement of each arm can be performed in 3 directions (x, y, z) plus rotation θ using a step by step motor. By taking the tank center as origin, the boundaries of each axis is given as follows :

           -1390 mm ≤ x ≤ 1390 mm -490 mm ≤ y ≤ 490 mm 0 mm ≤ z ≤ 900 mm 0 ≤ θ ≤ 360 • (2.40)
The minimum distance between 2 transducers that we can reach with such system is 100 mm. The theoretical accuracy of arm displacement and (x, y) directions is 100 µm that enables very accurate sound speed estimation.

Signal processing

For the signal generation and acquisition, we have used the National Instrument card PCI-6115 which features 12-bits multichannels analog to digital converter (ADC) and digital-to-analog converters (DAC) with sampling rate up to 2 MS/s. At the receiving side and before digitalization, signal coming from hydrophone is conditioned with a BK NEXUS charge amplifier. By using DAQmx drivers, the PCI-615 board is directly interfaced with MATLAB. The signal probe used for the experiment is a pure sinus wave centered at 150 kHz and sampled at 800 kHz with a duration of 6.6 µs. In order to estimate accurately the time delay of received pulse, the acquisition task is perfectly synchronized with a generated wave by using a digital trigger provided by the digital output pin of PCI-6115 board. The signal sent in the time and frequency domain is displayed in figure 2.15 and the recorded signal in the hydrophone is present in figure 2.16. The temperature is measured continuously to compute the sound speed in water.

Configuration

The configuration used is presented in figure 2.17. We form the horizontal antenna by moving the sensor 15 times by 1 cm. To have an accurate position of the source, as it is taken as a reference for the whole configuration, we measure the travel time for a round trip between the source and the plates to find finally h s = 16.18 cm. The distance between the source and the first hydrophone is 20 cm. The seafloor is simulated with two plates consisting of different materials. The first plate is made of Teflon, and the second is made of Resin. Their dimensions are 50 cm wide, 50 cm long and of different thickness (Table 2.6). The sound speed in the plates are found by using a contact measurement : knowing the thickness of the plates and using the travel time measurement between the boundaries of the plates we can estimate the sound speed. These measurements have been made using a transducer at 1 MHz center frequency [SL09]. The plates density is measured by simple weighting.

Parameters

Water

SCARAB 98 experiment signals

SCARAB is a measurement campaign that took place in June near the island of Elba in Italy (Figure 2.18). For more information, see to the reference [START_REF] Holland | High resolution geoacoustic inversion in shallow water : A joint time and frequency domain technique[END_REF]. The configuration is the configuration 1 presented in subsection 2.2.4. The source used is a 265 UNIBOOM of EGG model, an electro mechanical of 20 cm diameter, broadband (100 Hz -6 kHz) at 20 cm below the sea surface. The wave emitted by the source is reflected by the surface of water that inflected the directionality of the signal : this is the phenomenal of Llyod mirrors [START_REF] Jensen | Computational ocean acoustics[END_REF]. The grated part of the energy spectrum of the transmitted signal is below 3 kHz (50 cm wave length). The calculation of that with the point source model located 20 cm below the sea surface shows that their is no significant change in the directivity around the interest incidences angles (about 45 degree for the reflection on the first interface and less then 45 degree for the reflection on deeper interfaces). An example of the signal is displayed in figure 2.19.

Conclusion

In this chapter, simulated and real signals on which the effectiveness of the inversion methods will be tested are described. For simulated signals, the propagation equations of acoustic waves of a spherical source in an isotropic medium and an interface are presented. Then, the spherical waves are decomposed into a sum of plane waves with Sommerfeld integral. The reflection coefficient of plane waves on a fluid or solid multilayer has been studied. As a result, the simulated signals recorded on a layered seafloor composed of mitigating layers and a smooth interface are similar to the real signals. We have two types of real data :

-the first type is that of the SCARAB 98 experiment obtained with a vertical antenna and an impulse source at 200 m from the bottom -the second type consists of a tank experience obtained with a source and a receiver that moves to form an antenna in the case of impulse signal on two plates Teflon and Resin.

Chapter 3

The Image Source Method

Introduction

The aim of this PhD thesis is to characterize the seabed using acoustic methods. The characterization of the seabed by acoustic means lies in estimating the geoacoustic parameters from the signals measured by the instrumentation. With these instruments, the thickness, the sound speed and the attenuation are estimated. Their assessment is based on the inversion of a layered seafloor model. It is possible to have simpler and less costly geoacoustic inversion methods to classify the seabed. The image source method (ISM), a recently-developed method in S.Pinson's thesis [START_REF] Pinson | Sound speed profile characterization by the image source method[END_REF], is the method we seek to improve. In this chapter, we define the ISM which is used to find the sound speed and thickness of the layers.

Consequently, the first step is to present this method and how it is applied on a stratified medium. Next, how the positions of the image sources can be localized in constant sound speed space in order to get the parameters of arrival time and arrival angles. Then we present the algorithm for inverting the layer parameters (sound speed and thickness). We finish by presenting the model limits and error in the results of inversion on simulated data.

Representation of reflected signals with image sources

Principle

The theory of image sources is generally used in room acoustics [START_REF] Allen | Image method for efficiently simulating small-room acoustics[END_REF], for modeling the propagation in a wave guide [START_REF] Fawcett | A method of images for a penetrable acoustic waveguide[END_REF][START_REF] Fawcett | Complex-image approximations to the half-space acoustoelastic green's function[END_REF], or for modeling the reflection of radar antenna emissions on the floor [START_REF] Xu | An efficient analysis of vertical dipole antennas above a lossy half-space[END_REF]. Another application of the image sources is the validation of the performances of blind source separation algorithms, acoustic source localization, speech enhancement and channel equalization [START_REF] Lehmann | Reverberation-time prediction method for room impulse responses simulated with the image-source model[END_REF]. This method, emanating from geometrical acoustics, enables us to calculate the density of the sound energy in an enclosed or partially enclosed propagation medium. The image source model is intended to bring to a situation similar to that of a free field (which provides a very simple solution to the wave equation) by replacing each reflection on a wall by the corresponding image source. This implies, that each reflection is specular, and therefore that the wall can be considered perfectly flat and smooth with regard to the wavelength considered, so that the diffraction and the scattering can be neglected, and the curved walls can be considered as locally flat. In addition, it is common to equate the source to a point source assuming that the reception is at a sufficient distance, so that the waves emitted by the main source and each of the image sources are considered as spherical. In summary, the image sources models the waves reflected on the walls by a wave emitted by an image source symmetrical to the real source relative to the reflecting walls (Figure 3.1) which amplitude is affected by the reflection coefficient of the interface. The time for the direct and reflected waves to reach the reception corresponds to the distance between this point and respectively the main source or the corresponding image source in a two-dimensional space.

Water (c 0 ) Layer 1 (c 1 > c 0 ) x z Source (x s = 0, z s = 0) Receiver Image 0 x (0) I = 0, z (0) I = 2h S

Image sources for a layered medium

In this section, we present the application of the image source approach in the case of a spherical wave reflection by a layered medium under the Born approximation (Figure 2.1) instead of a wall or a waveguide. The main idea is that the reflection of an acoustic wave on a stratified fluid bottom (Figure 3.2) can be modeled with image sources. To do this there are some assumptions that must be verified to validate the approach of the image sources in the stratified medium. The ISM is applied to superficial sediments which are mainly composed of sand, clay and silt (section 1.2) which means they are composed of unconsolidated grains. The shear waves can be neglected because of their low sound speed and their high absorption coefficient in these environments.

Water (c 0 ) Layer 1 (c 1 , ρ 1 ) Layer 2 (c 2 , ρ 2 ) Basement (c 3 , ρ 3 ) x z Source (x s = 0, z s = 0) Receiver z K x K h s θ (0) 0 θ (0) 0 h r = h s -z K h s + h 1 θ (1) 1 θ (1) 0 θ (1) 0 h 1 h s + h 1 + h 2 θ (2) 2 θ (2) 0 θ (2) 0 h 2
The geological layers are parallel to each other : although the wavelength is chosen for a resolution of a few centimeters (λ = 15 cm for f = 10 kHz), the sedimentary structures range from centimeters to meters and have a low lateral variability in CHAPTER 3. THE IMAGE SOURCE METHOD compliance with the laws of Steno (section 1.2.2).

In our case, the influence of heterogeneities and the random roughness of the medium are not taken into account in this model because the signal reflected on a plane interface between two sediments will be slightly distorted since the frequency is low. But in the presence of strong heterogeneities and random roughness, the assumption of homogeneous layers is not valid.

2-Only first reflections are considered

Under the Born approximation [CS81, RF89], the amplitudes of multiple reflections are negligible compared to the first reflected one [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF]. The situation depicted in section 2.2.1 and sketched in figure 3.2 is the one we are interested in.

3-The incident angles on the interfaces are lower than the critical angle In order to remain in the frame corresponding to the previous assumption and to avoid phenomena such as head waves or total reflections on interfaces, a configuration is used in which the transmitter and the receiver are far enough from the bottom but are not too far from each other. Moreover, we should not be at normal incidence to have refraction. An incident angle of about 45 degrees is usually a good compromise.

Equivalent model

The bottom acts as a mirror showing a symmetrical image of the true source located on the same vertical axis as the source, and the same for each layers. Each reflection on each interface (Figure 3.3) can be identified by the receiver as an image source located in an equivalent medium : the structure (water + geological layer) above the interface and the reflection of this structure. Each image source is represented in a different system. The wave emitted by the image source passes through a layer whose thickness is doubled and the hydrophone is located at a height equivalent to the sum of the source height and that of the hydrophone from the bottom.

With this modelling, the reflections and transmissions in the layers are no longer respected so the values of the amplitudes no longer match the actual value. However, the arrival angles on the sensor and the travel time of the rays remain the same. As the order of the layers crossed has no effect on the travel time of the wave or on its arrival angle, we can rearrange the order and thus merge all the equivalent systems in a unique one where all layers thickness are doubled and where the image sources are localized at the interfaces (Figure . 3.4). Water (c 0 )

(c 1 , ρ 1 ) (c 2 , ρ 2 ) Basement (c 3 , ρ 3 ) x z Source Receiver 2h s h r = 2h s -z K 2h s + 2h 1 2h 1 2h s + 2h 1 + 2h 2

Geoacoustic inversion with image source method (ISM)

Real and virtual image sources

If the sound speed profile (SSP) is known, it can be taken into account and the image sources are exactly on the image interfaces and on the vertical axis below the source (black dots on Figure 3.5). But, in the context of geoacoustic inversion, this profile is what we are looking for and is thus unknown. So, instead of dealing with the "real" image sources, the recorded signals are modeled as being emitted by "virtual" images sources (VIS) which are the sources that give the good travel time and angle at an hydrophone location if all the medium is homogeneous like water, so when acoustic rays are straight lines. These VIS are no more located on the interfaces, neither on the vertical axis (+ symbols on Fig. 3.5 and 3.6). So, the virtual sources are detected by the hydrophone with an offset, which depends on the sound speed of the layers, from the vertical axis of the real source. But the arrival angle of the wave from the VIS on the receiver and the travel time of the wave between the receiver and the VIS remains the same as the real one.

Inversion principle (a) Parameters of ISM

For each layer i, we have a data pair (t (i) , θ (i) ) which are related to the geoacoustic parameters of the sediment structure. t (i) , θ (i) are the input parameters of ISM. The parameters to be estimated by ISM are the number of layers N , their thicknesses h i , and the sound speed of pressure waves c i .

Note : the travel time can be called arrival time when the recorded signals by hydrophones are synchronized with the emission time which is the case for all the signals processed in this thesis.

(b) The use of an array of hydrophones

From the signals recorded by hydrophone, we can obtain the travel time, but it is not possible to find the arrival angle by using a single hydrophone, we need at least two hydrophones. So the best way to obtain all the input parameters is to use an array of hydrophones in reception.

(c) Consequences of the use of an array

The locations of the real image sources (on the interfaces and the vertical axis below the source) are independent of the hydrophone position ; they only depend on the thicknesses of the seafloor layers (Figure 3.5). But this is not the case for the positions of the VIS : their locations depend on the hydrophone location. Therefore, when one uses an array to record the signals, one image source is not seen as a single dot but rather as a spot with a shape and a spatial extension that depends on the shape and extension of the hydrophone array (Figures 3.6 and 3.7). So, here one faces an inherent

3.3. GEOACOUSTIC INVERSION WITH IMAGE SOURCE METHOD (ISM) 57 Water (c 0 ) Layer 1 (c 1 > c 0 ) Layer 2 (c 2 > c 1 ) x z Source (x s = 0, z s = 0) Receiver z K x K h s Image 0 x (0) I = 0, z (0)
I = 2h S θ (0) 0 θ (0) 0 h r = h s -z K h s + h 1 Image 1 x (1) I = 0, z (1)
I = 2(h s + h 1 ) θ (1) 1 θ (1) 0 θ (1) 0 h 1 + Image 1 xI (1) , zI (1)
h s + h 1 + h 2

Image 2 x

(2)

I = 0, z (2)
I = 2(h s + h 1 + h 2 ) θ (2) 2 θ (2) 0 θ (2) 0 h 2 + Image 2 xI (2) , zI (2)
Figure 3.5 -Geometry for the "real" and "virtual" image sources in the case of an increasing SSP [START_REF] Guillon | Bayesian geoacoustic inversion with the image source method[END_REF].

problem of the ISM : to locate the image sources, one needs to use an array instead of a single receiver, but this implies that the image source is not a single point but a spot. The original situation sketched in figure 3.5 is in practice the one sketched in figure 3.6 : for an array of N h hydrophones, there are N h VIS instead of a single one for a given layer.

In order to define the data pair (t (i) , θ (i) ) for each layer, we need to define an equivalent receiver for the hydrophone array and also an equivalent VIS instead of dealing with the two spots. In the following, the median position of all the hydrophone positions is used to defined the equivalent receiver. This choice is made because the median is a statistic robust to outliers.

For the equivalent VIS, a new algorithm based on triangulation approach is used for this purpose (Chapter 5). It is not possible to find the theoretical position of the VIS in real data because we do not have the values of the SSP of the medium. However, in synthetic data and for a theoretical study of the model of ISM we can find the positions of the theoretical VIS. An example of the VIS in configuration 2 and Seabed 2 is presented in the figure 3.7. This model error is voluntary because it drastically reduces the number of input data to invert and therefore significantly increases the speed of inversion. A study of the model error and his effect in the inversion algorithm is presented in subsection 3.4.2.

Water (c 0 ) (c 1 , ρ 1 ) (c 2 , ρ 2 ) x z

GEOACOUSTIC INVERSION WITH IMAGE SOURCE METHOD (ISM)

(d) Algorithm

The ISM is based on the use of Snell-Descartes laws to invert the SSP of the seafloor from the arrival angle and the travel time of each equivalent VIS relative to the equivalent receiver. The ISM can then be summarized as a two-step algorithm (Figure .3.8) :

-step 1 : from the recorded signals, the number N s of virtual image sources is obtained and also the absolute travel time and arrival angle for each virtual image source i to the equivalent receiver. This step is presented in chapter 4. -step 2 : from these input parameters, the SSP is obtained by using the Snell-Descartes laws. This step is described in the next section 3.3.3.

(e) Advantage

The use of image source has a big advantage compared to the classical approaches (section 1.5.3) in term of computational time. Assuming that there are N -layers in the seafloor above an infinite basement, this will generate N s = N + 1 image sources. The reflected signals being recorded on an array composed of N h hydrophones, this leads to N t = N s × N h arrival times. This will be the number of inputs for a travel-time based inversion method (see e.g. [START_REF] Dosso | Uncertainty estimation in seismoacoustic reflection travel time inversion[END_REF]). Through the 1st point mentioned above, the ISM will convert this N t inputs into 2 N s inputs. This reduction in the number of inputs is the result of the direct model used here. It might imply a faster inversion algorithm than a global one but this might result in supplementary uncertainties that are investigated in section 3.4.

Geoacoustic parameters estimation

In this section, we assume that step 1 is already done and that we have already recovered the input parameters :

-the number of image sources N s -the travel times t (i) with i ∈ [0, .., N s -1] for each equivalent VIS to the equivalent receiver -the incident angle θ (i) i for each equivalent VIS at the equivalent receiver with respect to the vertical angle. As we know the input parameters, we are now able to calculate the sound speed and the thickness of each layer by a recursive method.

The first image source corresponds to the propagation of the signal only in water and enables us to determine the system geometry. Figure 3.9 presents the final equivalent system to model a stratified medium whose thickness and sound speed of layers are known to find the sound speed and thickness of layer l. 61 heq (0)=2h s -z K . Then we used the arrival time t (2) and the arrival angle θ

Water (c 0 ) Layer 1 (c 1 > c 0 ) Layer 2 (c 2 > c 1 ) Layer l (c l ) x z Receiver (x K ) S I O θ (l) 0 θ (l) 1 θ (l) 2 θ (l) l h r + h s 2h 1 2h 2 2h l
(2) 0 of the second VIS to find the parameters of the first layer. Hence the process is repeated for the following layers. So to generalize, knowing the travel time t (l) and the arrival angle θ (l) 0 , we can deduce the cl and heq (l).

To start, we compute the incidence angle in the layers by using θ (l) 0 to compute the incident angle from the first layer to the last one l -1.

θ (l) l-1 = sin -1 cl-1 sin(θ (l) l-2 ) cl-2 . (3.1)
Next, we compute the propagation time of the wave in the layer l

t ls = t (l) - l-1 n=0 heq (n) cn cos(θ (l) n ) . (3.2)
Then, we estimate the sound speed of the layer l

cl = c 0 sin(θ (l) 0 )t ls x K - l-1 n=0 heq (n) tan(θ (l) n ) , (3.3)
Finally, we estimate the thickness of the layer l hl = cl t ls cos arcsin cl sin(θ

(l) 0 ) c 0 . (3.4)
3.4 Limits and sensitivity of the ISM to bad input parameters detection.

To better understand the ISM, we treat in this part the limits and the different types of uncertainties that we may encounter in the different steps of the inversion. We start this section by a brief reminder of the uncertainties that were studied in S.Pinson's thesis [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF]. Then, two types of uncertainties are presented : the first is the one due to the model, the second is the one due to the bad detection of the position of the equivalent VIS which can cause an error in the input parameters.

Errors analysis [Pin11]

Like any other inverse method, the ISM has some uncertainties that may arise during implementation, hence generating uncertainties in the final results. In its thesis, the possible sources of errors were presented. There was also a study of their influence on the localization of the image sources and consequently on the SSP estimation. These errors belong to three different categories.

(a) Measurement errors

The positions of the source and the receivers are not recognizable, but this can be corrected by using the recorded signals and the emission time (see subsection 5.3).

CHAPTER 3. THE IMAGE SOURCE METHOD (b) Errors due to hypotheses

The geoacoustic hypotheses made for ISM (see subsection 3.2.2) can affect the results : the roughness of the geological interfaces can generate scattering waves which affect the detecting of arrival times and angles. The bathycelerimetric profile is generally constant in shallow water measurements. However, its variation can be a source of weak uncertainty.

The theoretical assumptions can also affect the results. For example, the multiple wave reflections between the interfaces are assumed too low in amplitude to be detected but if a multiple wave reflection is detected, this will imply the creation of a non existing layer in the estimated SSP. The wavefront is considered as spherical and the seafloor surface is considered as plane. However, the reflected wave is detected only in part by the antenna and the sphericity approximation of the wave in order to localize the image sources is certainly enough.

(c) Errors due to noise

The noise in the recorded signals can be the echoes due to multiple reflections between the interfaces, scattering phenomena, or noise in underwater environment. For the experiments on simulated signals with a noise ranging form 8 dB to -5 dB, a sensitivity of the estimated SSP to this noise has been noted. Some echoes can be submerged by very high noise. In this case the corresponding image sources are not detected and the layers will be mixed up.

(d) Conclusion

All these types of errors are still present in our case, except the error due to the mesh that no longer exists since we changed the first step of the detection algorithm (see section 4.2). However, another type of error related to the application of the model of image sources with a hydrophone antenna on a stratified medium is analyzed in this thesis. In the next section, we present this error which we is called the model error.

Study of model uncertainties

The principle of the ISM algorithm is based on the knowledge of absolute travel time and arrival angle between each equivalent VIS and a equivalent receiver. In order to test the extent of the uncertainty due to the choice of the equivalent VIS and hydrophone, we can use synthetic data because, in this case, it is possible to obtain the exact positions of all the VIS and thus, to study the influence of the choice of equivalent receiver and VIS.

We study the inversion results in the two geometric configurations 1 and 2 which are presented in chapter 2. But in this part we change the length of the antenna, in the configuration 1 to be 15 m with inter-distance of 1 m and its deepest hydrophone is located at z K =11.5 m from the bottom. This modification is done to not interfere the study of the parameters of length and number of hydrophone in the comparison. The To generalize the study, 3 different seabeds are used : the first seabed is the seabed 1, the second seabed is the seabed 3 and the third seabed is the seabed 4. For each configuration the geoacoustic parameters like thickness, the number of layers and sound speed are changed.

To begin, we calculate the theoretical time and angle between each image source and each hydrophone. The arrival angle of the first reflection is calculated as follows :

θ 1 K = tan -1 x K h s + h r (3.5)
with x K is the offset between the source and the K th hydrophone, h s and h r represent respectively the heights of the source and the hydrophone with respect to the seabed. To determine the incident angles in each of the remaining layers, we build a cost function that compares the horizontal distance x K with the same distance computed with the various incident angles :

(h s + h r ) tan(θ l K ) + 2 l=N l=2 h l tan arcsin c l+1 sin(θ l K ) c l -x K = 0 (3.6)
where h l is the thickness of layer l,c l is the sound speed of layer l, N is the number of layer and θ l K is the arrival angle of image source l to the hydrophone K. The optimization of this cost function is performed with MATLAB function "fsolve".

Finally the travel times are given by :

t l K = h s + h r c 1 cos(θ 1 K ) + 2h l c l cos arcsin c l+1 sin(θ l K ) c l (3.7)
Once the theoretical input parameters are found, the inversion is performed for both geometric configurations. The results are shown in tables 3.1 and 3.2. The first comment that we conclude from the inversion results in tables 3.1 and 3.2 is that the uncertainties found are very slight. The biggest uncertainty is around 2.99 10 -2 % in the sound speed and 7.53 10 -2 % in thickness in the two configurations which means less than 0.5 m/s and 3 mm in the worst case. So, the uncertainty of the model can be considered negligible for these configurations. But that does not prevent to make a comparative study between the two configurations.

Layers Error of

The second comment, we note that the thickness and the sound speed of layers can affect the uncertainties of the algorithm, but no general rule can be found there.

The third comment is that the uncertainties in configuration 1 are slightly above the configuration 2. This is probably due to the distance traveled by the wave : more than 150 m for the config.1 and 24 m for the config.2.

Finally, in configuration 1, we note that the uncertainties in the results, sound speed or thickness, becomes more important for the last layer than for the first one. We conclude that this algorithm becomes less efficient with depth which is not the case in the configuration 2.

Effect of length of the antenna or hydrophone number

To see the effect of length of the antenna or hydrophone number, the inversion is performed for the same geometric configuration 1 but with different antenna length. The same three different seabeds of the section 3.4.2 are used but with the configuration 1 of a vertical array made of 15 hydrophones with 64 m length and different inter-distances. The results are shown in table 3.3.

A comparison of the results reported in tables 3.1 and 3.3 shows that another parameter plays an important role in terms of the efficiency of the algorithm is the length of the antenna (or the number of hydrophones when we keep the same inter-distance and we add the number of hydrophones) used to find the position of the image source. So, we investigate the estimated parameters, SSP after finding the theoretical travel time and arrival angles, in the seabed 1, with 5, 10, 15 and 20 hydrophones. Figures 3.10 shows that the algorithm is more effective with 5 hydrophones than with 10 or 15. We move from 0.05 m/s and 0.09 m/s error in the first and the second layers with 5 hydrophones to 1.5 m/s and 2 m/s in the case of 15 hydrophones. The more antenna length is, the more virtual the image sources are extended, and the less accurate the equivalent image position becomes. To see more clearly the effect of antenna length, we plot the uncertainties values curve for different layers with different numbers of hydrophones for configuration 1 and seabed 1 (Figures 3.11 The figure 3.11 shows that the values of uncertainties are always greater for deeper layer even with a small number of hydrophones. It ranges from 0.05 m/s uncertainty in the first layer in the case of 5 hydrophones to 8 m/s at the last layer with 20 hydrophones. Nevertheless, the error values of the model do not exceed 8 m/s in the worst case which is quite acceptable.

LIMITS AND SENSITIVITY

To summarize, the uncertainties due to the model does not present an obstacle to the use of the algorithm and it is irrelevant to the type of antenna and the distance traveled by the wave. Even if the number of hydrophones remains an important parameter, its influence on the estimated results does not exceed the limits of uncertainty required.

Sensitivity of the ISM in the case of error in the image source position

The proposed algorithm is recursive. A measurement error on a layer can influence those below. To measure the extent of the sensitivity of this algorithm, we identify two possible types of errors. The first is an error in estimating the travel time and the second is an error in detecting the arrival angle.

In S. Pinson thesis [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF], a study of the sensitivity of the inversion algorithm to the 3.4. LIMITS AND SENSITIVITY OF THE ISM TO BAD INPUT PARAMETERS DETECTION. 67 position error detections is made : the localization error in distance and in the incidence angle. It has been concluded that the localization whether in distance or incidence angle propagates only to the next layer, i.e the errors on the next layer compensate those on the wrong layer and the errors in sound speed and thickness of the layers below are very low. In this part, this is reconsidered to see the influence of the choice of equivalent VIS and the effect of the length of antenna or the number hydrophones on the results.

Axial error

For axial error analysis the configuration 1 with seabed 1 is used. First, we calculate the time and the theoretical angle between each image source and each hydrophone as before. Then we add random values ∆t = 5 F s where F s is the sampling frequency, at each arrival time of the 3rd layer for each hydrophone. This step leads to advance or delay for the position of the image source (Figure 3.12).

Water (c 0 ) The equation 3.3 becomes as follows :

x z Equivalent receiver S SI 1 SI 2 SI 4 δ d SI 3
cl = c 0 sin(θ (l) 0 )(t ls + ∆t) x K - l-1 n=0 heq (n) tan(θ (l) n ) , (3.8)
The error on SSP one to an error in travel time is presented in figure 3.13.

Figure 3.13 shows the error of the SSP calculated for an error of travel time. This error causes a distance location error δd of the third source image in case the arrival angle is correct (Figure 3.12). We note that the error on the SSP focuses on the third and fourth layers. It is of the order ±15 m/s for layer 3 and 4, within ±1 m/s for the remaining layers. This shows that the error propagates only to the following layer by correcting the error of the previous layer's sound speed.

Radial error

As previously, for radial error we use the configuration 1 with seabed 1. We calculate the angles and the theoretical travel times. Next, we add a random perturbation value ∆θ of the order of 0.1 degree at each arrival time of the 3rd layer for each hydrophone without interfering with the arrival time (Figure 3.14). The equation 3.3 becomes as follows :

cl = c 0 sin((θ (l) 0 ) + ∆θ)t ls x K - l-1 n=0 heq (n) tan(θ (l) n ) , (3.9)
Figure 3.15 shows the error on SSP calculated for an arrival angle detection error of the third layer. This error causes a localization error in space δθ of the third image source in case the arrival time is correct (Figure 3.14). We note that the error in the SSP is concentrated only on the third and fourth layers. It is in the order of ±45 m/s for layer 3 and 4 and below ±1 m/s for the rest. This shows that the error spreads only to the next layer by correcting the error of the previous layer just like the error in the arrival time but the inversion algorithm is more sensitive to the error of the arrival angle.

Effect of number of hydrophone on the case of error on input parameters

As it has been seen before, the number of hydrophones used for detecting the position of image sources plays an important role in the uncertainty values of the algorithm. To complete this study, we compare the results of theoretical inversion with the inversions in different cases 5, 10, 15 and 20 hydrophones with axial and radial position errors in the third layer. .16 show that the algorithm is more effective in the case of 20 hydrophones than 10 or 5. We move from 20 m/s error in the third and fourth layers in the case of 20 or 15 hydrophones to 150 m/s in the case of 5 hydrophones. The more hydrophones we use, the less the error is small. Therefore, these results show that the number of hydrophones reduces the uncertainty on the inversion in the case of a detection error. Note that the error in the inversion parameters, SSP, is much more important and sensitive than the model error. To reduce the detection errors a large number of hydrophones must be used. Nevertheless, it should not be great to reduce the effect of the model error, 15 or 20 is a good compromise between the two.

Conclusion

In this chapter we have seen how to calculate the SSP of a fluid stratified environment with a description of the reflections on the seabed in terms of image sources. The information needed to locate are the travel time and arrival angles between each equivalent image source and equivalent hydrophone. The analysis of the signals received on the antenna enables us to locate the image source and thus to use the inversion algorithm based on Snell-Descartes law. We have also studied the model error of the method, the effect of the length of the antenna as well as the effect of the error on the input parameters, arrival angle and travel time on the inversed sound speed profile. The influences of each type of error and the number of hydrophones used have been studied

Chapter 4

Travel time detection and estimation

Introduction

As it has been shown in the previous chapter, the ISM is a geoacoustic inversion technique which provides the SSP of the sediment structure with a good resolution at a low computational cost. This method uses a broadband source and a vertical or horizontal array. The signal reflected by the seafloor is analyzed as being emitted by image sources relative to the real source with respect to the geological interfaces. The position of these images is determined by the SSP, their location is the key of the process. Using these positions as inputs, the SSP is obtained directly by using a simple algorithm based on Snell-Descartes refraction laws.

The present chapter deals with the first step of the algorithm of ISM, more precisely the detection of the arrival time. In the original work [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF], a migration map and a CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION semblance map are built from reflected signals and the combination of these two maps provides the image source localization to extract the input parameters. To avoid the construction of the two maps and their combination, a new strategy of low complexity based on the Teager-Kaiser energy operator is introduced. This operator amplifies the discontinuities and sudden amplitude changes in the signal while the soft transitions between samples are reduced. This property of the operator is used for the detection of the arrival time of image sources directly from the signal.

The chapter is divided into three parts : the first part presents briefly the original method to find the positions of the image sources. The second part is dedicated to the new method of travel time detection based on TKEO as well as all the solutions to improve the application. In the third part we present a combination of the TKEO and the time-frequency representation to improve the travel time detection in noisy signals.

Original ISM Algorithm for detection and localization

The original algorithm used for the detection and the localization of image sources is based on array processing and seismic techniques. It is thereafter labeled as ISM-MS. It was developed by Pinson [Pin11] and is described briefly in this section (See Appendix A for the details).

The idea is to migrate the recorded signals y K (t) of each the K th hydrophone in water to form a migration map I M . The signals focus on the image sources. Since the focus is not perfect, a semblance map I s is also computed which theoretically gives one at image source location and zero elsewhere. A mask is built on the semblance map where its value is above a given threshold and the maxima of the migration map are picked on this mask. This gives the localization of the image sources (Figure 4.1). Finally, the number of layers N s , the travel time t i and the arrival angles θ i are computed from these locations and the center of the array where i is the index of image source.

Migration map

The first step of ISM-MS consists in the reconstruction of the image sources by migrating the signals in a medium with constant sound speed (c 0 ). This migration focuses the recorded signals on the source and image sources positions. Accordingly, the signal received by the K th hydrophone is only the sum of the same pulses f (t) with different amplitude factors β m sent by all the monopolar sources M + 1. In the time domain, the signal is defined as follows :

ORIGINAL ISM ALGORITHM FOR DETECTION AND LOCALIZATION

75

S K (t) = f (t) * g(r s , r r K , t) + M -1 m=0 β m f (t) * g(r s m , r r K , t) + η K (t) (4.1)
where * is the convolution product, g(r s , r r K , t) is the Green function in the temporal domain, η K (t) is an additive noise assumed spatially white and f (t) is the transmitted signal. r s and r r K are respectively the position of the source and the receiver.

The back propagated signals are then filtered by a Gaussian window w(t) which has approximately the same duration as the signal f (t). Then the average over the N h hydrophones is calculated. The result (Figure 4.2), homogeneous to pressure, is mapped by :

I BW (r) = +∞ -∞ 1 N h N h K S bK (r, t)ω(t) 2 dt, (4.2)
where N h is the number of hydrophones. Since the image sources are located near the vertical source axis, a phase shift is applied on each recorded signal located at r r K to focus on the point r⊥(θ) located on this axis :

S H f K (d, θ) = D K (θ)s H K d c + τ K (θ) -τ c (θ) (4.3) where τ K (θ) = |r r K -r⊥(θ)| c 0 and τ c (θ) = |r r c -r⊥(θ)| c 0 present propagation times.
The exponent H is the Hilbert transform. The multiplication by D K (θ) = |r r K -r⊥(θ)| is the compensation of the spherical geometric divergence. d is the distance between the array center r r c and a point r of the image source map. Finally, the image sources are mapped with the following function :

I M (d, θ) = 1 N 2 h -N h N h n=1 N h q=1 S H f n (d, θ), S H * f q (d, θ) (4.4)
where the superscript * denotes the complex conjugate.

Semblance map

For easy automatic detection and a better localization of the image source, the semblance function is used. It consists of dividing I M (Equation 4.4) by the non-coherent summation of the focused signals :

I semb (d, θ) = 1 N 2 h -N h N h n=1 N h q=1 S H f n (d, θ)S H * f q (d, θ) 1 N h N h n=1 S H f n (d, θ)S H * f n (d, θ) (4.5)
This function takes the value 1 if all the focused signals are identical in phase and amplitude, and a value close to 0 if they are different from each other (Figure 4.3).

Interpretation of the ISM-MS

The main advantages of the ISM-MS are its low computational time compared to the classical inversion algorithms, good robustness against noise. In addition, the results obtained with real data are satisfactory for both vertical [START_REF] Pinson | Sound speed profile characterization by the image source method[END_REF] and horizontal antenna [START_REF] Pinson | Range dependent sediment sound speed profile measurements using the image source method[END_REF]. However, the ISM-MS is not completely automatic because the thre-4.3. NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 77 shold of the semblance map must be fixed beforehand to find the mask to be applied on the migration map. The automatic localization of the image sources in the migration map after applying the mask remains a delicate task. The ISM-MS algorithm is too complicated to be reproduced easily. The resolution of both semblance and migration maps depends on the mesh, so the finer the mesh the higher computing time. The computation time of the maps is very high compared to the global inversion algorithm. The false image sources can not be separated from the good ones. For these reasons, we try to improve this method by simplifying the algorithm of detection and localization of image sources in the homogeneous medium. The new approach is presented in the next section.

New detection algorithm based on the Teager-Kaiser energy operator

In this section, we present a new detection algorithm based on the Teager-Kaiser energy operator which is intended to estimate the arrival time of the waves reflected from the seabed directly on the signal in order to replace the map calculations used in the original method and to improve the resolution. After a review of the theoretical and physical foundations of this operator, we present the application of this operator on our simulated and real signals presented in chapter 2. We also present all the problems encountered in its implementation and their respective solutions proposed.

Time delay estimation problem

The problem of time delay estimation between noisy signals received at two spatially separated sensors has attracted considerable attention in various applications of sonar, radar, geophysics, acoustics, and biomedical engineering [START_REF] Carter | Coherence and time delay estimation : An applied tutorial for research, development, test and evaluation engineerings[END_REF]. Different algorithms with varying degrees of accuracy and computational complexity have been designed to address this problem. Many of these algorithms are related to the generalized cross-correlation approach introduced by Knapp and Carter [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF]. Some of the related schemes are Roth processor [START_REF] Roth | Effective measurements using digital signal analysis[END_REF], smoothed coherence transform [START_REF] Youn | On the roth and scot algorithms : Time-domain implementations[END_REF], phase transform [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF], and maximum likelihood estimator [START_REF] Carter | Coherence and time delay estimation : An applied tutorial for research, development, test and evaluation engineerings[END_REF].

Let x 1 (t) and x 2 (t) be two noisy signals received by two spatially separated sensors. To determine the time delay, the cross-correlation function is computed :

R x 1 x 2 (τ ) = E [x 1 (t)x 2 (t -τ )] , (4.6)
where E denotes the expectation operator. The argument τ that R x 1 x 2 (τ ) provides an estimate of the time delay. When the time shift matches the true delay correctly, the maximum peak value of R x 1 x 2 (τ ) takes the biggest value. Because of the finite time of x 1 (t) and x 2 (t), however, the cross-correlation can only be estimated. For ergodic processes, an estimate of R x 1 x 2 (τ ) is given by :

Rx 1 x 2 (τ ) = 1 T -τ T τ x 1 (t)x 2 (t -τ )dt , (4.7)
where T is the time observation interval.

CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION

Relation (4.6) relies on the stationary assumption of x 1 (t) and x 2 (t) and it is difficult to guarantee in practice. For non-stationary signals, cross-correlation function fails to give an estimate of the time delay between the two signals. The approximation given by equation 4.7 is valid provided that the integration time T is as long as possible without the signals become non-stationary. Thus, for short observation interval, cross-correlation function does not provide accurate time delay estimate.

To overcome these limitations, we investigate the feasibility of the Teager-Kaiser energy operator, which is not constrained by these assumptions, for time delay estimation. This quadratic energy operator is an instantaneous measure and thus it is independent of the time observation interval of the signals. Depending upon the instantaneous amplitude and instantaneous frequency of the signals, this operator is well suited for handling nonstationarity.

Definition of the signals energy

In signal processing, the energy of a signal is traditionally defined using the square of the signal magnitude, the envelope of squared signal magnitude, or the integral of squared signal magnitude. Time-frequency distributions are also used to determine the signal energy content in both time and frequency bins and particularly for the study of non-stationary signals. All these approaches estimate the signal energy without referring to what the signal represents [START_REF] Boudraa | Teager-kaiser energy operators in timefrequency analysis in time-frequency signal analysis and processing[END_REF]. An alternative approach is to define the energy of the signal based on physical principles. This approach is well evidenced in the speech signal production principle. Indeed, in nonlinear speech modeling, the dominance of modulation as a process in the production of speech points to the importance of analyzing signals from the point of view of the energy required to generate them [START_REF] Teager | Evidence for nonlinear sound production mechanisms in the vocal tract[END_REF]. It was also observed that in the conventional view of energy, tones at 10 Hz and 1000 Hz of the same amplitude have the same energy. However the energy required to generate the signal of 1000 Hz is much greater than that for the 10 Hz signal. One observation was that speech resonances can vary rapidly both in frequency and amplitude, essentially due to the separated air flow in the vocal tract [START_REF] Maragos | On amplitude and frequency demodulation using energy operators[END_REF]. This alternative notion of energy can be understood by focusing attention on the generation process of the signal [START_REF] Kaiser | On simple algorithm to calculate the 'energy' of signal[END_REF]. Based on the theory of oscillation of a simple spring-mass system, an energy operator is defined to measure the energy of a speech signal, which is produced by a nonlinear process. This concept of energy has further been extended to a large class of signals.

Teager-Kaiser Energy Operator

A simple form of the energy operator was derived by Kaiser, based on the findings of Teager and Teager [START_REF] Teager | Evidence for nonlinear sound production mechanisms in the vocal tract[END_REF], using basic physics of motion for a simple spring and mass oscillator [START_REF] Kaiser | On simple algorithm to calculate the 'energy' of signal[END_REF]. The obtained operator is called Teager-Kaiser energy operator (TKEO). Newton's law of motion of a mass m suspended by a spring with a constant K is described by a second order differential equation : K/m is the frequency of oscillation, and φ 0 is the initial arbitrary phase. The total energy E x of this undamped system in the Newton's physics is equal to the sum of the potential energy (in the spring) and the kinetic energy (of the mass) and is given by

d 2 x(t) dt 2 + K m x(t) = 0. ( 4
E x = K 2 x 2 (t) + m 2 d 2 x(t) dt 2 . (4.9)
By substituting d 2 x(t) dt 2 and x(t) in equation (4.9) we get,

E x = 2π 2 m(Af ) 2 ∝ A 2 f 2 , (4.10)
which is the true energy in the harmonic oscillator system. i.e., the energy required to generate the signal. We immediately see from equation (4.10) that the energy is proportional to both A and f . Note that according to equation (4.9), the conventional energy approximation x 2 (t) only takes into account the kinetic energy of the signal source.

Consider now a discrete-time signal representing the motion of the oscillator body :

x(n) = A cos 2πf f s n + φ o , (4.11)
where f is the analog frequency and f s is the sampling frequency. Using three equallyspaced samples x(n -1), x(n) and x(n -1), and appropriate trigonometric identities, it is easy to get

x 2 (n) -x(n -1).x(n + 1) = A 2 sin 2 2πf f s (4.12) Equation (4.12) is exact and is unique provided that the 2πf fs value is restricted to values less than π/2. This equation can be approximated if we notice that for small 2πf fs , sin( 2πf fs ) ≈ 2πf fs . The largest deviation from sin( 2πf fs ) is at f fs = 1/8. Thus, equation (4.12) can be reduced to

x 2 (n) -x(n -1)x(n + 1) ≈ 4π 2 A 2 f f s 2 (4.13)
Equation (4.13) holds provided that f s > 8f . The basic definition of the discrete TKEO is given by

Ψ[x(n)] = x 2 (n) -x(n -1).x(n + 1) (4.14)
A close look at equation (4.14) reveals that the TKEO involves nonlinear operation on the signal. Also, this operator is nearly instantaneous given that only three samples are required in the energy computation at each time instant.

The TKEO has the following properties [START_REF] Kaiser | On simple algorithm to calculate the 'energy' of signal[END_REF] :

-Ψ[.] is independent of the initial phase φ 0 , -Ψ[x(n)] = Ψ[x(-n)],
-Ψ[x(n)] is robust even when the signal passes through zero, as no division operations is required, -Ψ[x(n)] is capable of responding rapidly (in two sampling instants) to changes in both A and f .

CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION

In the continuous time domain the TKEO is defined as [START_REF] Kaiser | On teager's energy algorithm and its generalization to continuous signals[END_REF] :

Ψ[x(t)] = dx(t) dt 2 -x(t) d 2 x(

TKEO frequency Analysis

To understand the behavior of the TKEO, its output has been analyzed in the frequency domain [START_REF] Atlas | Quadratic detectors for general nonlinear analysis of speech[END_REF]. The TKEO belongs to the family of the quadratic operators defined by :

y(n) = L-1 l=0 x(n + l)x(n -l)h(l) (4.18)
where h(l) is an L-point set of quadratic filter coefficients. It is easy to see that the impulse response of the TKEO is given by

h(l) =      1 l = 0 -1 l = 1 0 otherwise
The frequency response of the filter is given by

H(e j2πf ) = e jπf 2(1 -cos(2πf )) (4.19)
Relation (4.19) indicates that the TKEO has the behavior of a high-pass filter and thus it accentuates the high-frequency content of the signal making it efficient, for example, for spike detection [START_REF] Mukhopadhyay | A new interpretation of nonlinear energy operator and its efficiency in spike detection[END_REF].

TKEO Application

As it has been shown, the TKEO is a nonlinear energy tracking operator and its output to a given signal, x(t), is the actual physical energy required to produce x(t) [START_REF] Boudraa | Teager-kaiser energy operators in timefrequency analysis in time-frequency signal analysis and processing[END_REF]. An important aspect of this quadratic operator is that it is nearly instantaneous given that, in its discrete version (Eq. 4.14), only three samples are required in the energy computation at each time instant. Equation (4.14) shows that TKEO computes a running estimate of the signal energy at each instant that takes into account the signal strengths at is immediate neighbors. Thus, it is justified that TKEO can be applied to the problem of detection [START_REF] Drira | Geoacoustic inversion with image source method : Improvement of the detection algorithm[END_REF] as it computes the instantaneous energy relative to its immediate neighbors. Also, it is evident from equation (4.14) the TKEO is expected to suppress slowly varying parts of the signal x(t) and to highlight abruptly changing parts (sudden bursts) of this signal. The TKEO has found many applications and particularly in speech processing [START_REF] Boudraa | Teager-kaiser energy operators in timefrequency analysis in time-frequency signal analysis and processing[END_REF]. Significant research on the theory and applications of the TKEO has been conducted [START_REF] Boudraa | Teager-kaiser energy operators in timefrequency analysis in time-frequency signal analysis and processing[END_REF]. The majority of the analysis in signal and image processing has mainly dealt with the properties of TKEO-based demodulation using the ESA and its extended versions and not with the operator itself [START_REF] Boudraa | Teager-kaiser energy operators in timefrequency analysis in time-frequency signal analysis and processing[END_REF]. In the present work we exploit the output of the TKEO for detection purpose.

The TKEO is directly applied on synthetic signals, of configuration 1 and seabed 2 with fluid layers, in two different cases : without noise (Figure 4.4) or with noise (Figure 4.5). For noiseless and moderately noisy signals the TKEO performs well. For very noisy signals the detection is not easy, thus the detection is improved by cross-correlation with the emitted signal (Figure 4.6). The detection process is explained in subsection4.3.7. TKEO. So the preprocessing improve the detection of travel time in noisy environment.

The operator is also tested on real data recorded in the tank experiment (section 2.3.1) with and without application of the cross-correlation. Since the recorded signals in the tank have high SNR, then the application of the crosscorrelation does not improve more the temporal signal resolution (Figure 4.7 (left)). But more false peaks are detected on the TKEO output of the recorded signal than on the pre-processed one (Figure 4.7 (right)) which makes the automatic detection of the arrival time difficult. For example in the second peak we detect 3 arrival time instead of one. These results show that even on less noisy data, the preprocessed signal has the same temporal resolution than the recorded signal, the preprocessing is a crucial step for the application of the TKEO to ensure an automatic and accurate detection.

The choice of the preprocessing

To estimate accurately the travel time it is preferable that the signal transmitted by the source is a zero phase which means that emitted signal is centered at t=0. In this 4.3. NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 83 case the travel times correspond to the maxima of the received echoes. When the emitted pulse is not zero phase, which is the case of our signals, the cross-correlation between the emitted signal and the registered signal can be used (Figure 4.8) :

s n (t) = s n (t ) * f * (-t ) (4.20)
where s n (t ) and f * (-t ) are respectively the signal recorded by the hydrophone and the conjugate of transmitted signal. The correlation operation extends the duration of the reference wavelet which decreases the temporal resolution of the signals. Consequently, it is preferable to use the method of Wiener deconvolution [START_REF] Sacchi | Statistical and transform methods in geophysical signal processing[END_REF] in order to conserve the temporal resolution.

s n (t) = T F -1 S n (ω)F * (ω) |F (ω)| 2 + W (ω) (4.21)
where T F -1 stands for the inverse Fourier transform. Sn(ω) and F (ω) are the Fourier transforms of s n (t) and f (t) respectively. W (ω) is the spectrum we wish to apply to the pulse and is a small constant value used to limit the impact of noise on the areas of the spectrum where there is less energy. In the following, the synthetic signals are processed with a matched filter since the temporal resolution of the emitted signal is well suited to the configuration of the medium and to the layers thickness. The real signals are processed with Wiener deconvolution since the propagation medium or the true shape of the signal sent are unknown.

Threshold

The TKEO-based detection is illustrated, as an example, on simulated signals of configuration 1 and seabed 1, corrupted with white Gaussian noise at SNR= 3 dB, and deconvoluted by a cross-correlation with the source signal. The peaks of the TKEO output (Figure 4.9) correspond to arrival times but can also be attributed to noise.

To overcome this problem, a threshold, T , is used to isolate the most prominent peaks that are considered as the reflection of the waves on the surface. More precisely, this threshold is introduced to minimize the missing of true peaks, while keeping the number of false detected peaks, caused by noise, within a reasonable limit. For decision, the following rule is used :

S ≥ T
The peak is taken into account,

S < T

The peak is not taken into account, (4.22)

where S are the amplitude values of detected peaks.

T value is taken as a scaled version of the average of the TKEO output :

T = C L L n=1 Ψ(n) (4.23) ψ(n) = 1 K K k=1 Ψ(x k (n)) (4.24)
where K is the number of hydrophones, L is the number of samples and C is a scale value which depends, essentially, on the SNR of the signal. There is no analytic solution for it but the following strategy is used to find C.

We recall that each prominent peak is associated to a given layer. For each SNR value and for each hydrophone, we identify the peaks corresponding to the arrival times (which are known for simulated signals) and we take the amplitude of the smallest peak, noted a K (m) where m is the SNR value. This operation is repeated for each hydrophone

NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 85

K and a function is constructed :

f (m) = min 1≤K≤N h [a K (m)]. (4.25)
Each point of the curve f (m) gives the smallest peak, corresponding to the interface between two consecutive layers, over all the K hydrophones. Note that if the impedance between layers is slowly varying, both peak detection and threshold selection are not easy. The curve in blue (Figure 4.10) corresponds to f (m) function. For decision, the following rule is used :

-If f (m) ≥ T , peaks corresponding to all layers are identified.

-If f (m) < T , some peaks cannot be detected and the associated layers cannot be identified.

For real data the SNR is unknown, C value is set to 1. If we take C > 1 there is a risk to detect peaks attributed to noise. If we take C < 1, very prominent peaks are well detected but some layers can be lost. For robust detection (against noise), we take the threshold that it is as close as possible to the curve of f (m). More precisely, the selected threshold avoids false peaks attributed to noise. As shown in figure 4.10, for m = 3 dB the best choice of C is 0.8.

To see the efficiency of the threshold value selection, tank data are used. Threshold values are calculated using equation (4.23) and the detected peaks are compared to true arrival times. This is possible, since for tank data both the layer number and the configuration are known.

Resolution improvement

The ISM is quite simple, but is still sensitive to errors of arrival angle estimation (subsection 5.2.2). Consequently, for better resolution of peak detection an over sampling of the deconvoluted signal is necessary (Equation 4.13) holds provided that f s > 8f ). But in practice oversampling of noisy signals can induce fluctuations on peaks which in turn complicates their detection (Figure 4.13). To remedy to this problem, a smoothing strategy based on polynomial filter is used.

Savitzky-Golay filters

Savitzky-Golay (SG) defined a class of filters well suited for smoothing and determination of derivatives at each point in a set of data equally spaced. These filters were designed to analyze the peaks in noisy spectra. The absorption peaks suggest which chemical elements are presented in the tested of materials [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF]. They showed that fitting a polynomial to a set of input samples and then evaluating the resulting polynomial at a single point within the approximation interval (window) is equivalent to discrete convolution with a fixed impulse response. It has been demonstrated, on noisy data obtained from chemical spectrum analysers, that these low pass filters reduce noise while maintaining the shape and height wave from peaks [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF]. let us consider a series of equally spaced noisy data points y(n). The SG filters work on the principle of local polynomial regression inside a sliding symmetric window [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF]. The idea is perform a polynomial approximation over the window, and assign the value of the concerned data point. This approximation is repeated for all samples of the series y. For each y(n), we least-squares fit a polynomial to samples y(-M ) to y(M ) inside the window, ε N , and calculate the coefficient of an N th order polynomial

ε N = M n=-M (p(n) -y(n)) 2 (4.27)
where a j are the coefficients of the polynomial function. The fitting optimized coefficients minimize ε N :

∂ε N ∂a l = M n=-M 2(y(n) -p(n))n l = 0 (4.28)
For better readability a matrix notation to describe this optimization problem is used. Thus (N + 1) equations of (4.28) can be stacked in matrix form as follows [START_REF] Schafer | What is a savitzky-golay filter[END_REF] :

a = (A T A) -1 A T y = Hy (4.29)
where y = [y(-M ), ..., y(M -1), y(M )] T and A = n i is the design matrix where n ∈ {-M, .., 0, .., M } and i ∈ {0, 1, ..., N }. Superscript T stands for matrix transpose. Invertibility of A T A is assured when the number of data points, 2M + 1, is greater than the order N of the polynomial. Thus we must have N ≤ 2M . Since just a 0 is required to obtain the smoothed version of y(n), the first row of H is multiplied by vector y to calculate a 0 . Once the window length is decided, the matrix A is known, therefore the first of row of H can be decided beforehand. Therefore, the coefficients of the filter can be calculated before the smoothing (filtering) is really taking place, which makes the filtering process very fast. The important observation is that the matrix H depends only on N and M and is independent of the input samples y [START_REF] Schafer | What is a savitzky-golay filter[END_REF].

An application of the SG filter to impulse synthetic signal is shown in the figure 4.14. The first point that can be noted from this study without interpolation and SG filtering is that the relative error averaged over 15 hydrophones on each layer does not exceed 1% in [0 dB -30 dB] interval (Figure 4.16). Beyond, the difference between the theoretical time and the calculated one becomes much more important. To ensure a visible interpretation of the figure, we stop at 0 dB. In the case of interpolation and SG filtering (Figure 4.17), the relative error, in the worst case, is of 0.01% up to -2 dB. In summary, these figures show that the relative errors are reduced when we improve the time resolution by the interpolation, which involves a very accrue detection of travel time.

The travel time is the key step in the new algorithm of ISM-TK by the fact that we use this parameter to find the position of the equivalent image source in the homogeneous medium in the next chapter by a very sensitive algorithm. So, a good detection of the travel time will move us forward to improve on the two input parameters of the geoacoustic inversion.

Improvement of matrix detection (a) Problem

To ensure the correct application of the method on real data at sea, our detection method is tested on SCARAB data signals (Figure 4.18). An additional problem appears : the different hydrophones do not detect the same number of peaks and some peaks detected by some hydrophones are not detected by the other. For example, hydrophone 1 detect 11 peaks whereas hydrophone 15 detects only 5 peaks.

NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 91

To solve this problem we propose a solution in the following section but before, we present briefly certain types of effects that can introduce this type of detection error.

(b) Possible causes

The detection can be affected by different physical problems. That the propagation of acoustic waves in the sedimentary medium can be affected in different ways, making the signal processing difficult.

When an acoustic wave strikes the interface separating two media with different acoustic impedances, a portion of the wave is transmitted while another portion is reflected by the interface. If the impedance contrast between the two layers is low, it can make it difficult for the various hydrophone signals to detect the peaks specially. This problem can be more important when the distance traveled by the wave is large because under the geometric attenuation the amplitude of the wave is lowered. So, if the distance between two different hydrophones is large, peaks can be detected by hydrophones and not by the last ones.

The arrival angles on the hydrophones are different, which affects the reflection coefficient directly related to the amplitude and the arrival angle. This makes the reflections visible on some hydrophones and invisible by others.

Directivity of the source may also make some peaks to be visible by some hydrophones and invisible for others Due to the large number of sound sources in the ocean, the underwater environment is noisy. The signal components considered as noise do not only come from ambient noise. In fact the echoes due to multiple reflections between the interfaces and the diffusion phenomena related to the roughness of the interfaces or inhomogeneities of the volumes will also involve a jamming signal. Some of them are so loud that they can be considered as a signal reflection on an interface and this can occur only for some hydrophones.

The shear waves "S" in the sediment medium are often ignored. In the study case, these waves parasite the signal because they originate from an energy loss at the level of the interface and create a delayed echo compared to the original response. The paper [START_REF] Guillon | Geoacoustic inversion of stratified seafloor with elastic layers by the image source method[END_REF] shows the effect of these waves on the ISM-MS and this describe and analyzed in section 6.2.

(c) Consequence of these problems on geoacoustic inversion

For a better comparison between the theoretical travel times and the detected ones, simulated signals of the configuration 1 and seabed 1 are used in order to known perfectly the theoretical times.

First the TKEO is applied on synthetic signals after preprocessing. The peaks detected are placed in matrix with N h = 15 lines and N s = 9 columns. Secondly, some errors are introduced into the travel time detection matrix of the echoes to simulate a poor detection under the effect of the problems cited above. Then, we make a comparison with the inversion values in case of all the detection values are taken into account. From this test, if we miss a peak in the matrix detection, it initially appears that the program replaces the values missing in the time detection matrix by the immediate neighboring values. The lack of a peak for a hydrophone hence shifts all the line values, which made a completely false matrix.

In the first case, some peaks have been removed on some hydrophones. For this, we taken a random example : hydrophone 1 has not detected the layers 4 and 7, hydrophone 2 has not detected the layers 6, hydrophone 4 has not detected the layer 2 and hydrophone 7 has not detected the layer 5.

The sound speed errors for each layer show the influence of a detection failure (Figure 4.19). In this case, the detection failure affects mainly the layers 7, 8 and 9. An error of about 15 m/s has appeared on the layers 1 to 6 which proves that the algorithm remains efficient in the case of small detection case, but it is much more bigger on layer 7 which is reduced a bit on layer 8 but doubles over the last layer. The error of sound speed in layer 9 is around of 300 m/s but it is not illustrate in the figure because we have zoomed to see the other errors.

In the second case, the third peak of the signal is withdrawn on consecutive hydrophones (2, 4, 8, 11, 12 and 14). This error at layer 3 results in an error of about 15 m/s (Figure 4.20). As we have seen in the chapter 3, the error in a layer propagates only to the next layer, so the errors on the layer 4 compensate those on the 3 layer. This error is then tiny for the following layers (layers 5, 6, 7 and 8) till it is accumulated in the last layer. It is important to note that a peak corresponds to the echo reflected by a layer. Detection errors for some peaks or for a layer on some hydrophones cause sound speed and thickness errors at the level of inversion. It is therefore imperative to correct such a defect created by the absence of a value by developing an automatic detection method of this defect.

CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION (d) Correction Principle

To solve the problem of different numbers of peaks detected by the hydrophones, an alignment to the first peak is performed for all the TKEO output. This step allows to remove the difference in travel time between the hydrophones for the first image source. However, this step is not sufficient because we still are not able to identify the peaks that corresponds to each layer as illustrated in figure 4.21. The difference is more noticeable for example for the 10 th and the 15 th hydrophones. The same phenomena is observed on TKEO output of Scarab data (Figure 4.22). For efficient for peaks detection (reduction of the false negatives), the first peak that has the shortest travel time is taken as a reference, regardless of the hydrophone. Then, the signals of all hydrophones are checked to determine whether there is a peak detected in an time interval computed for each hydrophone. This time interval is related to the order of the hydrophone and that of the layer ; indeed, a larger time gap is expected between the first and the fifteenth hydrophone than between the fifth and the tenth hydrophone. A matrix is created to identify whether the peaks are detected or not the example of first case is shown in the table 4.1. The matrix elements are 1 or 0. The peak detection is noted by 1 and absence of peak by 0. We turn the counter back if the hydrophones have not detected the peak to avoid going beyond the peaks. The correction principle of the matrix detection is finally incorporated into a loop to process the entire signal.

CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION

Hydrophone / Layer 1 2 3 4 5 6 7 8 9

1 1 1 1 0 1 1 0 1 1 2 1 1 1 1 1 0 1 1 1 3 1 1 1 1 1 1 1 1 1 4 1 0 1 1 1 1 1 1 1 5 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 1 1 1 1 7 1 1 1 1 0 1 1 1 1 8 1 1 1 1 1 1 1 1 1 9 1 1 1 1 1 1 1 1 1 10 1 1 1 1 1 1 1 1 1 11 1 1 1 1 1 1 1 1 1 12 1 1 1 1 1 1 1 1 1 13 1 1 1 1 1 1 1 1 1 14 1 1 1 1 1 1 1 1 1 15 1 1 1 1 1 1 1 1 1

Proposed solutions to the inversion

Different corrections are possible to fill the missing values in the detection matrix :

-A-Insert the neighboring value ; -B-Replace the missing value with the average value of neighboring values ; -C-Take into account only the hydrophones where the peak has been detected. The results of each type of correction are presented in figure4.24 for test case.

For correction A, an oscillation of the speed error which is not reduced is observed. This error oscillates between ±5 and ±20 m/s over layer from 2 to 9.

For correction B, the error is limited between -20 m/s and 10 m/s over layer 2 -9 with a strong tendency to remain in the order of 5m/s. Finally, for correction C, gives almost zero error over layer to 3, of the order of 2 to 3 m/s for layer 4-7 and of the order 8 m/s for layer 8 and 9.

The correction giving the lower errors is the one that eliminates the hydrophones where the peaks have not been detected at the inversion. But this is less efficient when the number of hydrophones "defective" becomes too high.

TKEO associate with time-frequency representations

The performances of the inversion method can be affected in very noisy environments. Due to the sensitivity of the TKEO to high level noise, spurious peaks can be produced which provide wrong time delay estimates. This is also the case in the presence of impulsive noise as encountered in underwater acoustics. A solution to this problem can be the denoising of the received signals on the hydrophones. However, the denoising probably makes the expected denoising effects unsatisfying due to over-denoising or under-denoising. The over denoising can affect the temporal structure of arrivals and under denoising does not reduce the number of false peaks. For more robustness to noise, we propose to perform the detection of peaks in the time-scale or in the time-frequency domains. This is essentially motivated by the fact that sub-band signals derived, for example, from a wavelet decomposition are narrow band where the TKEO is more efficient for peaks detection than on the original signal which is can be wideband.

Time-frequency representations

It has recently been shown the interest of time-frequency representations (TFRs) for peak detection or source separation [START_REF] Yilmaz | Blind separation of speech mixtures via timefrequency masks[END_REF]. In this work, we investigate the interest to combine the TKEO and the resolution of a TFR such as the spectrogram or the scalogram for peaks detection and particularly in noisy environments. One advantage of these two TFRs is that they are free of cross terms. The Stockwell transform which is CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION closely related to the spectrogram and to the scalogram is also used. For the illustration of the time-frequency approaches, the data used in this chapter is the synthetic signal (Figure 2.10) computed with configuration 2 and seabed 2 (subsection 2.2.4). An example of synthetic signal corrupted with an additive white Gaussian noise with SNR set to -9 dB is shown in figure 4.25.

Short-Term Fourier Transform (STFT)

The Short-Term Fourier Transform (STFT) is a TFR widely used for studying non-stationary signals. It is obtained by applying the Fourier transform by a fixed-sized, moving window to input signal. The STFT gives a 2D spectrum for time-frequency analysis and is defined by

Y s (t, f ) = +∞ -∞ y(τ )h(τ -t)e -j2πf t dτ, (4.30)
where h(τ -t) is a sliding window of a fixed length and y(τ ) is the analysed signal. The energy density surface of the STFT, called spectrogram, is given by

S p (t, f ) = |Y s (t, f )| 2 . (4.31)
The STFT has a fixed resolution in time and frequency, in the same manner as the Fourier transform. The frequency resolution ∆f = 1/N ∆t is controlled by the length of the window h(t) where N is the number of samples. The spectrograms of the signal with or without noise are displayed in figure 4.26.

Continuous Wavelet Transform (CWT)

As the STFT, the continuous wavelet transform (CWT) is a well established tool for time-frequency signal analysis. It allows a better accurate localization in time of high frequency signal features compared to the STFT [START_REF] Addison | Time-frequency analysis of biosignals[END_REF]. The CWT has gained attractiveness in representing and preserving the signal energy in the presence of noise. The CWT of a signal y(t) is defined as :

T (a, b) = 1 √ a +∞ -∞ y(t)Φ * t -b a dt, (4.32)
where Φ * (t) is the complex conjugate of the analyzing wavelet function Φ(t). T (a, b) represents the wavelet coefficient where a and b are the dilation and location parameters respectively. The contribution to the energy of y(t) at the specific a scale and b location is given by the 2-D wavelet energy density function known as the scalogram (Figure 4.27) :

E w (a, b) = |T (a, b)| 2 . (4.33)
The relative contribution to the total energy of y(t) at (a, b) is given by the variance of the wavelet defined for the CWT as the scale-dependent energy distribution by

E(a) = 1 C g +∞ -∞ E w (a, b)db, (4.34)
where C g is the admissibility constant. In practice E(a) is used to determine the dominant energetic scales within y(t) (Figures 4.28

Stockwell Transform (ST)

The continuous S-transform (or Stockwell transform) of a signal y(t) is defined as [START_REF] Stockwell | Localization of the complex spectrum : the s transform[END_REF] : The signal y(t) is windowed (or multiplied point by point) with a scalable window Gaussian function where the standard deviation is proportional to the period :

S T (τ, f ) = +∞ -∞ y(t) |f | √ 2π e - (τ -t) 2 f 2 2 e -i2πf
σ(f ) = 1 |f | (4.36)
As the time domain width of the window is determined by the frequency f , it can easily be seen that the window is wider in the time domain for the lower frequencies and narrower for the higher frequencies. In other words, the window provides good localization in the frequency domain for the low frequencies, whereas it provides good localization in time domain for higher frequencies.

The ST has a time localizing Gaussian window that is translated and the Fourier kernel (different from wavelet kernel) that is stationary which localizes the real and the imaginary components of the spectrum independently. The Fourier kernel selects the frequency being localized and hence it localizes both the amplitude and phase spectrum. Thus, it retains an absolute phase of the signal which is not provided by the CWT.

The ST is closely related to the CWT and the STFT. More precisely, the ST has a form similar to the STFT except that the width of the window varies with the frequency, and can be considered as the CWT with phase-correction. Besides, the ST used time-frequency axis rather than the time-scale axis used in the CWT. Therefore, the interpretation on the frequency information with the ST is more straightforward than with the CWT which is beneficial to remove high frequency noises. The linear property of the ST ensures that for the case of additive noise, the ST gives :

TKEO in time-frequency domain

TKEO-STFT (TK-STFT)

The spectrogram is considered as an image in which the peaks are detected. The application of the TKEO on the spectrogram can be viewed as an enhancement operation of this image representation and thus improves the robustness of the peaks detection. More precisely, the TKEO is applied along the rows (time) and columns (frequency) of S p (t, f ) and the filtering results summed as follows :

Ψ[S p (t, f )] = Ψ t [S p (t, f )] + Ψ f [S p (t, f )] (4.38)
where Ψ l [.] is the operator applied along to the l th component. The created time-frequency binary mask for time delay estimates is calculated as follows :

M p (t, f ) := 1 if Ψ[S p (t, f )] ≥ th 0 otherwise (4.39)
where th = mean[mean(Ψ[S p (t, f )])] is an automatic threshold used to isolate the most prominent peaks that are considered as the reflections [START_REF] Drira | Image source detection for geoacoustic inversion by the teager-kaiser energy operator[END_REF]. The associated scale-dependent energy distribution, E(a), is illustrated by figure 4.28. The peak of E(a) is located at a = 8 which corresponds to the dominant energetic scale. Thus the peak detection is performed on this scale. The TKEO is applied to component of scale a = 8 and the result is illustrated in figure 4.31. This figure shows that the TK-CWT reduces the number of false detected peaks compared to the TKEO method where the operator is applied directly on the y(t) signal.

TKEO-ST (TK-ST)

The TFR of ST is seen as an image in which the peaks are detected. As for the TK-STFT, the application of TKEO on this image improves the robustness of the peaks detection [START_REF] Drira | Estimation des temps d'arrivée pour la méthode des sources images en combinant l'opérateur de teager-kaiser et la transformée de stockwell[END_REF]. Specifically, the operator is applied along the lines (time vectors) and columns (frequency vector) of S T (τ, f ) and the filtering results is summed as follows :

Ψ[S T (t, f )] = Ψ t [S T (t, f )] + Ψ f [S T (t, f )] (4.40)
As in TK-STFT the binary mask is given as follows :

M T (t, f ) := 1 if Ψ[S T (t, f )] ≥ ths 0 otherwise (4.41)
where ths = mean[mean(Ψ[S T (t, f )])] is an automatic threshold used to isolate the most prominent peaks that are considered as the reflections [START_REF] Drira | Image source detection for geoacoustic inversion by the teager-kaiser energy operator[END_REF].

Figure 4.32 shows the binary mask M T (t, f ) extracted from the ST representation of y(t) with SNR=-4 dB. The peaks are well evidenced on the mask with no false peaks.

Results

Three methods are compared, namely, TK-CWT, TK-Spectrogram and TK-ST on synthetic data with configurations 2 and seabed 2. The theoretical travel times are compared to those detected estimated with the three approaches on synthetic data for different SNR values and for 15 hydrophones. Then, the relative error in % is calculated for each approach in different cases : depend on the scales and the better estimates are obtained for the dominant energetic scales (a ≥ 8) and this confirms the findings of the figure 4.28.

CHAPTER 4. TRAVEL TIME DETECTION AND ESTIMATION

Some remarks

Some remarks can be pointed out concerning the combination of the TKEO and the TFRs (CWT, Spectrogram and ST) for the time delay estimation in noisy environment. From the results obtained on synthetic signals the following conclusions can be drawn : TK-CWT is effective to detect the desired peaks, even for low SNR (-9 dB), and it reduces the number of false alarms. The computational load of the TK-CWT is lower than that of the two other methods. However, the results of the TK-CWT are conditioned by the choice of the adequate pre-determined mother wavelet and the scale parameter. TK-Spectrogram and TK-ST are efficient and perform similarly for SNR above -4 dB but the TK-ST requires slightly more computational time than the TK-Spectrogram. The main advantage of the TK-ST over the TK-CWT and the TK-Spectrogram, is that it does not require any input parameter. Compared to the TK-CWT, the TK-Spectrogram requires only one parameter, the size of the analyzing sliding window. Finally, for the three methods, the threshold for the binary mask extraction is computed automatically from the data.

Conclusion

To improve the performance of the ISM and to reduce its computational load a new approach of travel time detection based on the TKEO is presented. After recalling the basics of the original ISM, the principle of the TKEO is presented. The properties of this quadratic operator are highlighted showing their interest for peaks detection. The effectiveness of the TKEO has been shown on synthetic and real data. To identify the peaks corresponding to the arrival times from the output of the TKEO a binary decision strategy has been developed. For noisy signals, the TKEO based detection method does not perform well. Thus, a smoothing of the signals using the Savitzy-Golay filter is used. This filter is well suited for noise reduction while maintaining the shape and height of the noisy signal peaks, which are usually washed by a moving average. Under the effect of some disturbances, such as seabed noise, shear waves or differences in the spatial positions of hydrophones, some peaks detected by a hydrophone cannot be detected by another. To overcome this problem, an algorithm for the correction of the detection matrix has been developed. Due to the sensitivity of the TKEO to high-level noise, spurious peaks can be produced, which provide wrong time delay estimates. As a solution to this problem, the detection of the peaks is done in time-scale or in time-frequency domains where the TKEO is applied to the coefficients at each scale or in each line and each column of spectrogram presentation. The detection results obtained for very noisy signals show the interest to combine the effect of the TKEO with a time-frequency (time-scale) representation as a strategy for peaks detection in noisy environment.

Chapter 5

Image sources localization and inversion results

Contents

Introduction

As it has been seen in chapter 3, the ISM is a geoacoustic inversion technique that is based on the positions of the equivalent image sources. Their location is the key of the process to find the input parameters which are the travel time and the arrival angle. Using these parameters as input, the SSP is obtained directly by using a simple algorithm based on Snell-Descartes refraction laws.

In the last chapter, we have described the process to obtain the first parameter which is the travel time of each image source to each hydrophone of the antenna. In this chapter, we focus on the localization of the equivalent VIS in a homogeneous medium. The positions of this equivalent VIS gives the input parameters which are the arrival angle θ (i) and the travel time t (i) . Then, by using the inversion algorithms the seabed parameters are found.

The chapter is divided in three parts : the first part deals with localizing the image source by using a triangulation approach. The sensitivity of the algorithm in the case of error is analyzed and a solution is proposed to correct it. Then in the second part, two methods to correct the true shape of the antenna are presented as well as their application on synthetic and real data. Finally, we present the results of the inversion algorithm on the synthetic and the real signals presented in chapter 2.

Image sources localization for input parameters estimation

As it has been seen in chapter 3, the real image sources are "seen" by the hydrophone array as a spot of virtual image sources (VIS). The input parameters of ISM are, for each layer, the travel time t (i) and the arrival angle θ (i) between the equivalent VIS and the equivalent receiver which is defined as the median of the hydrophone positions. The next step of ISM is so to define the equivalent VIS. To do so, we propose to use the travel times between the VIS and the hydrophones estimated in chapter 4.

To better explain the principle let consider one particular layer i. We have the t (i) K travel times corresponding to the propagation from the virtual image sources of this layer to the hydrophones. These travel times are converted in distances

D (i) K with K ∈ [1, N h ]
by the knowledge of water sound-speed c 0 where N h is the number of hydrophones in the antenna. For one hydrophone K, its corresponding virtual image source is located on a circle with a radius equal to D K . But the position of the VIS in this circle is unknown. The proposed solution is to use the intersections of circles to localize an equivalent image source close to the equivalent image source found by the median of VIS. This strategy is called triangulation and is explained in detail in the next subsection.

Triangulation principle

Let start by a simple intersection of two circles. In this case, we can find the positions of the point M(x M , z M ) by searching for the intersection of the circles (x K i , z K i , R i ) and (x K j , z K j , R j ), where (x K i ,z K i ) and (x K j ,z K j ) represent respectively the sensor coordinates in the (x, z) plane for the first circle i and the second circle j. (R i , R j ) are the radii of the circles which represent the distances between the point M and their respective hydrophones.

The intersection of two circles is a former maths problem, based on the resolution of the following system of equations :

(x M -x K i ) 2 + (z M -z K i ) 2 = R 2 i (x M -x K j ) 2 + (z M -z K j ) 2 = R 2 j (5.1)
In order to simplify the calculation, we suppose that the two hydrophones are vertically aligned (x K i = x K j ) (figure 5.1). After solving the problem, the coordinates of the point M are found as follows :

z M = R 2 j -R 2 i + z 2 K i -z 2 K j 2(z 2 K i -z 2 K j ) x M = R 2 i -(z M -z K i ) 2 .
(5.2)

IMAGE SOURCES LOCALIZATION FOR INPUT PARAMETERS ESTIMATION

109 Two solutions can be found for these equations. To find the right solution, two options are possible : either use a third hydrophone, either keep only the solution in the part of the space we are interested in. Considering the physics of the problem considered here, we can choose the second option.

Water (c 0 )

x z A 2 A 1 A 3 M R j R i x M z M x M z M Figure 5
.1 -The triangulation to find the coordinates of the point M

Analysis of the sensitivity

In this part, we make a mathematical study of the uncertainty of the position of the image source in case of error in the travel time detection or error in the hydrophone positions. This is done by adding an uncertainty ε i on the radii as follows :

R i = R i + ε i R j = R j + ε j (5.3)
We suppose that the errors are small compared to the distance traveled by the wave ε i,j R i,j . Thus neglecting the second order terms (ε 2 i = 0) leads to :

R 2 j -R 2 i = R 2 j -R 2 i + 2(ε j R j -ε i R i ).
(5.4)

The coordinates of the point M can then be written as follows :

x M = R 2 i -(z M -z K i ) 2 + 2R i ε i z M = R 2 j -R 2 i + z 2 K i -z 2 Hj 2(z 2 K i -z 2 K j ) = z M + ε z , (5.5)
with

ε z = ε j R j -ε i R i z K i -z K j .
(5.6)

We note that the ε z error is especially related to the radii of the circles R in the numerator and to the inter-sensor in the denominator. So, the distances traveled by the waves and the inter-sensors have a big effect on the sensitivity in the case of error in the travel times. The bigger the radius and the smaller the inter-sensor, the larger the errors.

This result shows a very high sensitivity to errors because ε z is not negligible. As an example, if the distance R i is to order 10 2 and i is to order 10 -1 , which means that the error is 0.1% in the radius, the uncertainty ε z becomes 10 (10 % error on z M ) which is not acceptable.

In the next subsection, we propose a solution to reduce this sensitivity to be able to use the triangulation algorithm for the image source detection.

Solution

To ameliorate the solution, we compute all the intersections of the N h × (N h -1) circles, pair by pair, to find all the coordinates (x i,j M , y i,j M ) possible. This is done by solving the non-linear system :

(x K i -x) 2 + (z K i -z) 2 = D 2 k i (x K j -x) 2 + (z K j -z) 2 = D 2 k j (5.7) where K i ∈ [1, N h ] and K j ∈ [1, N h ] K i = K j .
This system is solved by the fsolve Matlab function with Gauss-Newton algorithm. Obviously, this system has two solutions. We retain the solution which lies near the position of the real image source, i.e. near the nadir of the real source. There is also a cleaning step in the algorithm to get rid of the unrealistic solutions. Finally, the position of the equivalent VIS is obtained by computing the median of all the intersection solutions previously found (Figure 5.2). In the end, to apply the inversion algorithm, the travel time t (i) and the arrival angle θ (i) for layer i are computed between the positions of the equivalent VIS and the equivalent receiver.

To see the efficiency of the equivalent image source selection with the median of all the circle intersections, the geometric configuration 1 and seabed 1 are used. The choice of this configuration is made to test the proposed solution on a simulated configuration in which the time traveled by the wave is much larger than the inter hydrophones distances. The coordinate values of the three first image sources are located using the triangulation algorithm, by taking the median of all the intersections, in the case of different error orders on the radii. The detected positions are compared to the median of the equivalent VIS. The results are presented in table 5.1. Table 5.1 -The positions of image sources in the cases of different order of add errors ε compared to the theoretical equivalent VIS.

From these results, it appears that the sensitivity of the triangulation algorithm is much lower with the median of all the intersections. If we take the same example as before, the distance R i is to order 10 2 (150 m in this configuration) and is to order 10 -1 (10 cm), the uncertainty ε z is 0.22 m in the worst case which is acceptable. But if is to order 1 m which is the case of bad hydrophone positions, the uncertainty ε z become 1.8 m. To see the effect of this radius errors we present the inversion results in figure 5.3.

The errors on the sound speed profile are to the order of 60 m/s in some layers which are not acceptable. Hence, we must reduce the effect of this high sensitivity of the triangulation method to be able to use it. So, firstly we increase the accuracy of travel time detection by oversampling the signals (section 4.3.8) and secondly we correct the uncertainties of hydrophone positions. The correction approach is explained in the next section.

Deformation and correction of the array

In general, the geoacoustic inversion methods make assumptions about the propagation environment and about the geometry configuration. As one of these inversion methods, the ISM is very sensitive to the acquisition geometry, particularly the geometry of the receiving antenna. Basically, the array of hydrophones used is often assumed linear throughout the experimentation. But this is not necessarily checked in the case of a long, flexible antenna and towed in a seismic flute type and in the case of an experiment in the sea such as the one used to support this study. Because of the sea currents, the horizontal or vertical antenna may be deformed. Similarly, a dip in the structure may imply that the antenna is not perpendicular or parallel to the seabed.

Various estimation methods related to the distortion errors of the antenna have been proposed in the literature [START_REF] Marcos | Calibration of a distored towed array using a propagation operator[END_REF][START_REF] Edwards | Numerical techniques for efficient sonar bearing and range searching in the near field using genetic algorithms[END_REF]. For example in [START_REF] Marcos | Calibration of a distored towed array using a propagation operator[END_REF], two reference sources, with known relative positions are used to reduce the distortion and correct the positions of the sensors. In [START_REF] Bourennane | Locaoperator of wideband sources with estimation of an antenna shape[END_REF] the authors assume that the shape of the antenna is known. In [START_REF] Pinson | Range dependent sediment sound speed profile measurements unsing the image source method[END_REF], the correction of the antenna is based on the comparison of the pulse phases from the source and the first reflected wave on the bottom in order to determine the arrival time of the error.

In this subsection, we focus on the estimation of the shape of the array. We propose new methods to estimate this shape based on the TKEO which is used to estimate, with great accuracy, the travel times of the direct path and the first reflected path [START_REF] Drira | Estimation de la forme d'une antenne d'hydrophone à l'aide de l'opérateur d'énergie de teager-kaiser[END_REF].

As in [START_REF] Branderbug | A point mechanical model for the dynamics of towed arrays[END_REF][START_REF] Edwards | Numerical techniques for efficient sonar bearing and range searching in the near field using genetic algorithms[END_REF], the correction is reduced to an optimization problem.

To estimate the true shape of the antenna, we detect firstly the arrival time of the direct path and the first reflected path recorded by the sensor network. This processing step is based on the quadratic energy operator TKEO, as we have seen before in chapter 4.

Secondly, the distance between the source and the hydrophones (Fig. 5.4) are calculated simply by :

D s K = t s K c 0 , (5.8)
D im K = t im K c 0 , (5.9)
where c 0 is the sound speed in water which is considered constant, t s K and t im K are respectively the travel times of the direct path and the first reflection on the bottom to the hydrophone K. We propose two algorithms, explained in the next sections, to detect the real shape of the antenna with the advantage of fast calculations and easy implementation. The first one consists in solving the nonlinear system of the triangulation problem. The second one consists in optimizing a function of the squared errors.

Water (c 0 ) Interface x z Source (x s = 0, z s = 0) K 1 z K1 x K1 K 2 z K2 x K2 K 3 z K3 x K3 h s Image 0 x (0) I = 0, z (0)
I = 2h S D im 1 h r2 h r3 h r1 D im

Triangulation

If the rectilinear antenna undergoes small deformations, it is possible to apply corrections to the hydrophones coordinates in order to deduce the actual geometry of the antenna. This geometric reconstruction of the shape of the antenna is possible from the knowledge of the source coordinate and the inter-hydrophones distances of the antenna. This is done by solving the non-linear system (equation 5.10) to find the parameters xK , zK , and θ(1) 0K .

(x K -x K-1 ) 2 + (z K -z K-1 ) 2 = D 2 int 2h s -zK c 0 cos( θ(1) 0K ) = t 0k (x K -x (0)
I ) 2 + (z K -z (0) I ) 2 = xK (2h s -zK ) tan( θ(1) 0K ) = (D im K ) 2 (x K -x s ) 2 + (z K -z s ) 2 = (D s K ) 2
(5.10) where h s is the height from the source to the bottom, θ(1) 0K is the arrival angle of the first reflection in the bottom and D int is the distance between two adjacent hydrophones. This system is solved by the fsolve Matlab function with Gauss-Newton algorithm in an iterative way. At the end of this step, we have the three parameters (x K , zK , and θ(1) 0K ) which are the hydrophones coordinates as well as the arrival angle of the first image source that presents a main parameter to estimate the geoacoustic parameters.

Optimisation error approach

The second correction approach is based on the optimization of an error function. For this, the distance found in equations 5.8 and 5.9 between the sources and the hydrophones are used, rewritten here as :

D s,im K = (x s,im -x K ) 2 + (z s,im -z K ) 2
(5.11)

where the index (s) and (im) are respectively the source and the first image source. The idea is to add error terms δ xK and δ zK on the original positions of hydrophones. Then, we seek to minimise there errors. In the presence of errors, we have :

x K = x K + δ xK ,
(5.12)

z K = z K + δ zk , (5.13)
with δ xK x K and δ zk z K . Thus, the distance becomes :

D s,im K = (x s,im -x K -δ xK ) 2 + (z s,im -z k -δ zK ) 2 , (5.14)
which, after a limited development of the first order becomes :

D s,im K = a 2 s,im + b 2 s,im 1 + a s,im δ xK + b s,im δ zK a 2 s,im + b 2 s,im
(5.15) with a s,im = x s,im -x K and b s,im = z s,im -z K .

δ xK and δ zK are estimated by optimizing the cost function (equation 5.16) :

argmin (δ xK ,δ zK ) N h K=1 D s,im K -a 2 s,im + b 2 s,im 1 + a s,im δ xK + b s,im δ zK a 2 s,im + b 2 s,im
.

(5.16)

We find the following results :

δ xK = D s,im K a 2 s,im + b 2 s,im a s,im - a 2 s,im + b 2 s,im a s,im - b s,im a s,im δ zK .
(5.17)

δ zK = D s,im K a 2 s,im + b 2 s,im b s,im - b 2 s,im + a 2 s,im b s,im - a s,im b s,im δ xK .
(5.18)

As the following equations are bijective, the correction can be found as one following the other. So, if the correction following δ xK is known we can deduce δ zK or vice versa.

To determine the error on one of the axes we use the data of the direct path as well as the reflected one. The following relation is found :

δ zK = D im K √ a 2 im +b 2 im b im - D s K √ a 2 s +b 2 s bs + √ a 2 s +b 2 s bs - √ a 2 im +b 2 im b im b im a im -bs as (5.19)
Finally, the δ xK is deduced with the equation 5.17.

Results

(a) Applied correction algorithms in synthetic data

To test the efficiency of the proposed algorithms, we have used the configuration 2 (horizontal antenna see subsection 2.2.4). A set of data is simulated with a distorted antenna : to simulate a deformation in the shape of the antenna, the hydrophones positions are randomly distributed by adding random values in the original position (x K , z K ) in order of 1 m each time. Then, the travel times between the source and the first image source to each hydrophone are computed. Finally, we apply both correction algorithms to find the true shape of the antenna. For the triangulation approach, we do not need a priori information about the positions of the hydrophones except the inter-sensor that is supposed constant. For the optimisation error approach we have started by assuming that the antenna is straight. The results of the correction are presented in figure 5.5.

For both approaches the results show that the shape of the antenna is found very close to that of the true shape. Exceptionally, the theoretical positions are perfectly found by the triangulation method, so in figure 5.5 the red dots are indistinguishable and the uncertainty does not exceed 10 -7 m in the worst case. For the optimized approach the uncertainty is very small for the first hydrophone positions in order of 10 -4 m but for the last one the uncertainty is in the order of 10 -2 m.

To test the efficiency of these approaches in more realistic conditions, we have added a perturbation in the order of some samples to each travel time detected. Then, the two correction methods are applied as before to find the antenna shape. The results are presented in figure 5.6. The shape of the antenna is found close to that of the true shape but less accurate than before. The uncertainty is in the order of 10 -2 m for both of them. The triangulation approach is more sensitive than the optimisation error approach in the case of travel time error.

To conclude, both correction methods estimate correctly the true shape of the antenna in the case of random perturbation in the space and with an error associated Figure 5.5 -Shape of the antenna corrected with the triangulation and optimisation error approaches. Blue : the horizontal antenna, Black + : antenna shape found by the optimisation error approach, Green : antenna shape found by triangulation approach. Red : theoretical position.

Figure 5.6 -Shape of the antenna corrected with the triangulation and optimized methods in the case of a perturbation in the order of three samples to each travel time detected. Blue : the horizontal antenna, Black + : antenna shape found by optimized approach, Green : antenna shape found by triangulation approach. Red : theoretical position.

to the travel time detection. Each approach has advantages compared to the other. The triangulation approach is more accurate to find the antenna shape but it is more sensitive to the error associated to the travel time detection. However, the triangulation approach is more practical by the fact that we do not need a priori information about the positions of the hydrophones.

(b) Effect of true shape correction on inversion results

To see the effect of the antenna correction on the inversion algorithm to estimate the SSP of the seabed, we present the results of the inversion algorithm on synthetic and real signals. For the synthetic data the optimisation error approach is used to show the efficiency of this algorithm on a horizontal antenna and essentially because we have a priori information about the antenna shape. For the real SCARAB 98 signals the correction is done by the triangulation approach because we do not have a priori information about the antenna shape.

Synthetic data

The acquisition geometry is the configuration 2 (see subsection 2.2.4) and the seafloor is seabed 2. The geoacoustic parameters are estimated before and after the correction. The positions of the hydrophones are perturbed by adding random values between 0 and 1 m. We note that the inversion of the data is strongly improved by the knowledge of the geometry of the antenna. Before the correction of the antenna shape, the errors are around 10 % in the sound speed and 16 % in the thickness. But after the correction of the antenna shape, the geoacoustic parameters are correctly estimated.

Real data

For the SCARAB data, in S.Pinson thesis [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF], a simple correction of the antenna inclination relative to the horizontal under the effect of the slope was done. But this is not sufficient in our algorithm because the image sources are very poorly detected with the assumed linear antenna. This is probably due to the sensitivity of the triangulation algorithm. For this reason, after travel time detection, we apply the triangulation approach to find the true shape of the antenna without priori information. An example of the antenna shape found for real data (SCARAB 98 data) is shown in figure 5.7.

First, we note that the distance between the corrected antenna and the initial one is around 19 m, which explains the poorly detected image sources before the correction. Second, we note that the shape of the corrected antenna is inclined at the top around 2 m.

For the inversion results, in the case where the antenna is assumed to be rectilinear at 206 m above the source, the results are absurd with unrealistic values. That is why they are not presented. The inversion results obtained, after the correction of the antenna shape, are compared to those of Holland and Osler [START_REF] Holland | High resolution geoacoustic inversion in shallow water : A joint time and frequency domain technique[END_REF], which have already proven their reliability for the first 6 meters with direct measurements on sediment care (Table 5.3). As reported in the table, correcting the shape of the antenna improves the quality of the inversion results in real data. The inversion results are close to Holland and Osler results, which proves that despite the sensitivity of the triangulation algorithm, the correction of the antenna shape allows us to use the ISM-TK. The inversion results confirm the entire algorithm of ISM-TK namely the travel time detection with TKEO, the triangulation approach and the approaches of the correction of the antenna shape. In the next section we present the inversion results of ISM-TK on the signals presented in chapter 2. These inversion results are also compared to the results of the original approach ISM-MS.

Inversion results

In this section, we present the results of the ISM-TK algorithm to estimate the SSP. The data inverted are first the synthetic ones, and then the real ones.

Synthetic data

Synthetic data with a horizontal antenna

We start with seabed 2 and configuration 2 (see subsection 2.2.4). The estimated geoacoustic parameters, from the signals by the ISM-TK and ISM-MS methods, are presented in table 5.4. In this example, 15 hydrophones have been used to estimate the geoacoustic parameters with ISM-MS. For the ISM-TK, firstly 15 hydrophones have been used in order to have the same conditions as ISM-MS. Secondly 2 hydrophones have been used because to have a better resolution and less uncertainties in the results of synthetic data, as we have seen in chapter 3, a lower number of hydrophones is used which means a short length of antenna.

Parameters

Input These results show that, overall, the proposed method ISM-TK performs better in terms of estimation of thickness and sound speed parameters compared to the ISM-MS method. This is mainly the case when the length of antenna is short. Furthermore, the computation time is reduced by a factor ten in this configuration.

The advantage of the ISM-TK algorithm is that it is fully automatic because the threshold is automatically computed compared to the ISM-MS which need to set the threshold by the user in the semblance map to find the mask. Moreover, the ISM-TK is more straightforward to implement compared to ISM-MS. In ISM-MS, the resolution depends on the mesh. This means that to have good results, the mesh needs to be small enough which implies a larger computation time. This is not the case for ISM-TK where there is no need to set any parameters.

In the case of perfect knowledge of the geometric configuration, with no uncertainty in the hydrophone positions, only 2 hydrophones are used by ISM-TK to localize the image source and to have perfect inversion results. In ISM-MS the antenna-processing needs more hydrophones to have a good localization of the image sources.

Synthetic data with a vertical antenna

The second data set tested come from seabed 1 and configuration 1 (see subsection 2.2.4). The estimated geoacoustic parameters, from the signals by the ISM-TK and ISM-MS methods, are presented in table 5.5. Firstly, in this example 15 hydrophones are used to inverse the geoacoustic parameters in order to have the same conditions as Scarab 98. Secondly, 5 hydrophones are used to inverse the geoacoustic parameters with less uncertainty.

Parameters

Input The main points that can be deduced from this table are : -The uncertainties in the thickness are very slight in both methods, less than some millimeters. -The ISM-MS has not detected the layer number 6. So, the layers 5 and 6 are combined. The thickness is the sum of the thickness of the two layers and the value of the sound speed is an average of the sound speed of these two layers. -The ISM-TK detects all the layers. With ISM-TK a few uncertainties are noted in the case of 15 hydrophones and with no uncertainties in the case of 5 hydrophones. -ISM-TK reduces the computation time by a factor fifty in this configuration compared to the ISM-MS. To sum up, the results in the table show that, as in the previous case, the ISM-TK performs better in terms of estimating the thickness and the sound speed parameters compared to the ISM-MS, especially when the length of antenna is short. At the level of the computing time, the ISM-TK algorithm is much faster and more accurate than ISM-MS mainly if we increase the number of layers e.g moving from factor ten to factor fifty if the number of layers in the seabed change from 2 to 9.

Tank experiments

To investigate the feasibility of the proposed method and compare the performance of different time delays detection on real data, experiments have been carried out in water tank. The principle is to reproduce in a controlled environment a configuration similar to what one might have in the sea : a point source, a hydrophone, an array and a layered seafloor (section.2.3.1). The advantage is to benefit from an ideal laminated floor whose characteristics are known. The IR of the system composed of plates is then obtained after a Wiener deconvolution as pre-processing. The signal obtained after the Wiener deconvolution is shown in figure 5.8. In addition, the temporal resolution is good enough to visualize three reflections on the panel system interfaces (arrivals 2, 3, and 4).

ISM-TK algorithm with tank data

The ISM algorithm is applied on the signals recorded in the tank with a horizontal array of sensor and provides the thickness and sound-speed for the two plates. The SSP is obtained with the ISM-TK by using the TKEO to find the travel time. Then, the triangulation approach is used to find the equivalent image sources to extract the input parameters which are the travel time t i and arrival angle θ i of each layer i. A correction of the antenna shape is done by the triangulation correcting approach. Finally, the inversion results are presented in table 5.6.

The results are compared with the results of ground truth. The comparison is very satisfactory. In particular, the sound speeds are very close to the ones measured directly at 1MHz : the difference is below 0.4% in the first layer and 0.5 % in the second. The estimation of the thickness is less accurate for the first layer is around 30% but still satisfactory for the second less than 0.8%. The problem of thickness estimation of the first layer should be address in a future work.

Method

TK-TF algorithm with tank data

The recorded signals in the sensor of the tank have a good time resolution. The SNR of the signals is not measured but visibly is very high. As it have been seen in the last subsection the use of TKEO is enough to find accurate travel time. But to check the validity of the new approach that combined the TKEO and Time-frequency approaches for application to noisy signals, the travel time of the tank signals are detected with the different TK-TF approaches. The SSP obtained with the ISM using the new approach (TK-CWT), (TK-spectrogram) and (TK-ST) are presented in Table 5.7. The results are compared with the results found with the TKEO applied directly on the signals and to the ground truth. The comparison between the different approaches and the ground truth are satisfactory for some of them and not for the other. In particular, the use of TK-CWT have slightly improved the inversion results of the sound speed compared to TKEO but this is not the case in the thickness. The inversion results of TK-Spectrogram and TK-ST are close to the ground truth but less accurate than TKEO.

In summary, the idea of combined the TKEO and Time-Frequency approaches is to increase the capacity of the algorithm to detected the travel time in very noisy signals which is not the case of the tank signals. The application of these approaches on real data shows that the travel time detected with TKEO directly on the signals are better than the combination approaches except the TK-CWT that have the same resolution in 5.5. CONCLUSION 123 real data. As a conclusion of the TK-TF, the best approach in noisy or in clean signals and for both synthetic or real signals is TK-CWT.

Inverse results of real data Scarab 98

Estimation results of the ISM-TK, the ISM-MS and Holland & Osler method [START_REF] Holland | High resolution geoacoustic inversion in shallow water : A joint time and frequency domain technique[END_REF] are compared and are reported in Table 5.8. Since the ground truth, provided by coring, is only available for the first fifteen meters of the sediment, comparisons are restricted to this depth. Since Holland & Osler method compares favorably against the ground truth [HO00], it is used as the reference. Overall, our results confirm the findings of Holland & Osler (Table 5.8). However, only four layers are identified by ISM-TK method. Also, the first layer has not been identified. This can be explained by the complexity of the sedimentary stratification of the first meters. The main advantage of the use of ISM-TK is the speed of the results computing. We can find the inverse results in real time, it is 50 times more quickly than the ISM-MS which is considered faster than the other inverse algorithm.

The recorded signals in the hydrophone of the SACARB 98 data have a medium time resolution. As it have been seen in this subsection the use of TKEO is not enough to find accurate the travel time of all the reflection waves on the layers. It is interesting to test the new approach that combined the TKEO and Time-frequency approaches in this real data. This is what will be done in a future work.

Conclusion

In this chapter, we have seen how we use the travel time detected in chapter 4 to find the second input parameter which is the arrival angle by the triangulation approach. Based on triangulation technique, the high sensitivity of the algorithm is studied. So, to validate this approach on real data, the position errors of hydrophones are well corrected by two different methods before the use of triangulation step. Finally, we have presented the results of ISM-TK inversion in both synthetic and real data which are compared to

Introduction

The ISM provides the SSP using a two-step strategy. First, the signals reflected by each layer are detected for each hydrophone and the arrival time of signals are estimated by the TKEO. Then, using the triangulation method after a correction of the hydrophone position, the arrival angle between the image sources to each hydrophone on the antenna are estimated. In the second step, the SSP is computed via the Snell-Descartes laws of refraction using the two input parameters, travel time and arrival angle, found before. In this chapter, we extend the ISM for the estimation of other geoacoustic parameters namely the sound speed of shear waves, the density, and the absorption coefficient of the sediment layers.

The chapter is divided into three parts : the first part presents the approach used to estimate the sound speed of shear waves. The objective of the second part is to present CHAPTER 6. EXTENSION OF ISM TO OTHER GEOACOUSTIC PARAMETERS ESTIMATION two methods to find the density and the absorption of the sediment layers. In the final part, to check the validity of these approaches, the results of inversion are compared to the theoretical values for simulated data and to the ground truth of real experiments in a tank which is performed to obtain data in controlled conditions.

6.2 Shear wave sound speed inversion

Shear wave in stratified medium

In the ISM, the seafloor is modeled as a layered fluid medium, with flat, horizontal, and homogeneous interfaces. The estimated parameters are the number of layers, their thickness, and the sound speed of compressional acoustic waves inside each layer. As it has been seen in section 1.4.1, the assumption that the sediments are fluid media is very often made in the geoacoustic community and it is reasonable in many cases, especially for soft superficial sediments. But, rigorously, marine sediments can support shear wave (see e.g. [START_REF] Jackson | High-frequency seafloor acoustics[END_REF] and references therein).

As the shear wave is not included in the ISM model and therefore in the inversion algorithm, its existence may affect the obtained results. In [START_REF] Guillon | Geoacoustic inversion of stratified seafloor with elastic layers by the image source method[END_REF], we have concluded that the presence of shear waves inside the seafloor affects the accuracy of the ISM with errors on the sound speed of compressional wave and uncertainty on the layer's thickness. This is supposed to be an effect of mode conversion inside the sediment stack. But, image sources due to the shear wave are clearly seen and with an appropriate modelling and algorithm, shear proprieties may be inverted. For this reason, the aim of this section is to quantify the shear waves inside the sediment stack on the ISM performances.

We start our study with geometric configuration 2 (section 2.2.4) which is tested with seabed 2 parameters (table 2.3) : the first geoacoustic configuration uses only one fluid layer between the water and the basement (by simply putting to 0 the shear sound speed) and the second configuration uses one solid layer between the water and the basement.

The Figure 6.1 shows the propagation of acoustic waves in the case of a one elastic layer over a fluid basement. Due to mode conversion, the arrival of a P-wave on a solid layer may generate 4 types of waves as follows : the wave PPPP (trajectory number 1), the wave PSPP (2) that arrives at the same time at the receiver with the wave PPSP (3), the wave PSSP (4). The wave with the (M) index is the multiple reflection.

In figure 6.2, the signal corresponding to the configuration of a fluid medium (indicated in blue) shows three peaks : the first is the reflection on the first layer, the second is the reflection on the bottom and the third one which is the lowest corresponds to a multiple reflection. For the signal corresponding to the solid configuration (indicated in red), we also notice three peaks : the first is the reflection on the first layer, the second is the reflection on the bottom which is lower than in the case of a fluid medium because the energy of the wave is divided into P and S waves and the third is due to the shear wave. In this case, the multiple reflection is negligible since the P wave is already quite low from the first reflection. (1) (M )

T 12 T 23 T 12 R 23 R 21 T 23 h 1 θ 1 θ 1 θ 2 θ 1 θ 2 θ 3 (

CHAPTER 6. EXTENSION OF ISM TO OTHER GEOACOUSTIC PARAMETERS ESTIMATION

In the case of two solid-layers the problem of mode conversion is more complicated. The signal simulated in solid medium is presented in figure 2.11. The arrival of a P wave on two solid layers can give 16 types of waves. Some of these waves arrive at the same time as others. It is still very difficult to detect these waves because of their low amplitudes and mainly as we can confuse them with multiple reflections which have the same orders of magnitude in this case. Consequently, the analysis presented in this section is restricted to the case of a single elastic layer.

Detection of shear wave with TKEO

The TKEO is applied on the signals recorded to find the arrival times corresponding to the reflections of P longitudinal waves and S transverse waves, at the interfaces (Figure 6.3). We note that the arrival times of the first and second reflections are the same as in the case of a fluid medium, the third peak (in blue) is the wave (PPSP and PSPP) and the third peak (in green) corresponds to the multiple reflection for a fluid layer. This can be seen more clearly when the image sources are drawn in the homogeneous medium using the triangulation approach (Figure 6.4).

In the case of a solid layer, the positions of SI 1 and SI 2 are the same as the fluid layer. Since the S shear wave is slower than the P wave and because the position of image source in a homogeneous medium is related to the sound speed, the image source corresponding to the shear wave is located behind the vertical axis of the real sources. This is observed for the image source SI 3 S. SI 3 M is on the same vertical line in the case of a fluid layer since they have the same sound speed in the layer.

As a consequence, the location of the image source relatively to the vertical axis allows us to distinguish the shear waves from the multiple reflections. This spatial separation of image sources presents a advantage of the ISM in order to estimate the presence or absence of shear waves in some sediments.

Estimation of shear wave sound speed

After detecting the positions of image sources, SI 1 and SI 2 , we can extract the travel time and the arrival angle of each one. Then, the thickness and the P wave speed of the first layer are inverted as shown in Chapter 3.

To calculate the sound speed c s of the S wave, we draw the equivalent diagram in figure 6.5 which presents the equivalent system to model a propagation of acoustic wave with conversion mode in the first layer whose thickness and compression sound speed c p are known.

Since the thickness, the sound speed c p and the arrival angle of SI 3 are known, the distances x t1 and x tp in the figure 6.5 are calculated as follows : x tp = h1 tan(θ

x t1 = (hs + hr) tan(θ (S) 0 ), ( 6
(S) P ), (6.2)
where h1 is the estimated thickness of layer 1 and θ (S) P of the P wave in the first layer is :

θ (S) P = arcsin c1 c 0 sin(θ (S) 0 ) . (6.3)
These two distances provide the distance x ts :

x ts = x h -(x tp + x t1 ). (6.4)
The angle of the S wave in the first layer is obtained by : θ (6.6)

The results of inversion are shown in the table 6.1 :

The inversion result is close to the theoretical value, the uncertainty is around 1%. The uncertainty is probably caused by two effect : To estimate the density, we use the local reflection and transmission coefficients. For a plane wave arriving in a medium i at a flat interface spearing this medium to medium j, these coefficients are written :

R ij = Z j -Z i Z j -Z i and T ij = 2Z j Z j -Z i = 1 + R ij (6.7)
where

Z n = ρ n c n cos(θ n )
is the normal impedance for this wave in medium n.

Estimation of reflection coefficients

Reflected amplitude

At first, it is necessary to find the corresponding amplitudes of the waves reflected by the layers from the recorded signal without any pre-processing. For this purpose, we use the TKEO, but as we have seen before in section 4.3.6 the resolution of TKEO output is very poor without pre-processing. So, we use the arrival time detected with the pre-processed signals. Then, we search the amplitude max in the recorded signal in the area of arrival time detected before by TKEO, in a time interval equivalent to the duration of the wave sent. These processes are repeated for each reflection.

Amplitude of the emitted signal

To find the value of the reflection coefficient, the amplitude A s of the emitted signal is needed. In our configuration the sound speed is supposed constant and the positions of the sensors are reliable. In this case the amplitude B 0K of the direct path to hydrophone K is easily and accurately measured. To calculate A s , we need to estimate the losses (geometrical and attenuation in water) which are given by the following relation :

dL p = 20 log 10 (D 0K ) + α(D 0K ) (6.8)
where D 0K is the distance between the source and the K th hydrophone and α is the absorption coefficient. In our case, due to the low frequency in real experience Scarab and the small distances in the tank, we assume α=0.

Hence the amplitude of the emitted signal is obtained from the amplitude B 0K as follows :

A s = B 0K 10 dLp 20 , (6.9)

The reflection coefficient generally depends on the impedances and absorption coefficient of the layer on which the wave is reflected. To determine the density of the first layer which corresponds to the first reflection coefficient, we study the effect of the absorption coefficient in the reflection coefficient.

Figure 6.8 shows that the absorption coefficient has a negligible effect on the values of reflection coefficients within the range of used angles (less than the critical angle). But this may be not the case for the other layers with the embedding inside the sediment since the refraction angle increases.

Estimation of the density of the first layer

In this paragraph, we are interested in finding the density of the first layer. For this the parameters available are the amplitude of the source A s , the amplitude of the CHAPTER 6. EXTENSION OF ISM TO OTHER GEOACOUSTIC PARAMETERS ESTIMATION reflected wave B 1 , the sound speed c1 , the thickness h1 and the arrival angle θ 1 0 . The first step is to calculate the distance D 1 traveled by the wave between the source and the equivalent receiver :

D 1 = h s cos(θ (1) 0 ) + h r cos(θ (1) 0 )
(6.10)

where h s and h r are respectively the heights of the source and the receiver and θ

(1)

1 is the angle in the first layer which is found by the Snell-Descartes law. After that, the reflection coefficient is estimated by the following relation :

R12 = B 1 D 1 A s = Z 2 -Z 1 Z 2 + Z 1 (6.11)
where Z 1 and Z 2 are the impedance of water and of the first layer. The impedance of the first layer is then given by :

Z 2 = Z 1 1 + R12 1 -R12 . (6.12)
Finally, the density of the first layer is estimated by :

ρ2 = Z 2 cos(θ (1) 1 ) c 2 (6.13)
The results of this estimation are present in the section 6.4.

The estimation of the densities of the remaining layers is based on the calculation of refection coefficient of stratified medium. This calculation uses two different strategies :

-By following the individual trajectory of the wave in the stratified medium (subsection 6.3.4), -By considering the overall layered medium (subsection 6.3.5). These two strategies lead to two different methods for density estimation.

The amplitude method (AM)

In this section, we present a method based on the wave propagation and the local reflection coefficient to find the densities of the remaining buried layers [DMB + 15] with an iterative way. The parameters available are the amplitude of the source A s , the amplitude of the reflected wave B iK , the sound speed ci , the thickness hi and the arrival angle θ i 0K where i is the index of the layer number. As we have seen before, we assume that the absorption of the first layer does not have an effect on the local reflection coefficient. As for the layers below, we firstly estimate the densities assuming that the medium is without absorption in order to simplify the calculation. This is the case of weak absorption in a homogeneous medium. Secondly, we estimate the densities taking into account the absorption of layers in order to have more accurate estimation of the density and to have access to the absorption coefficient.

Case without absorption

In this paragraph, we present how we can find the density of layer l in the case without absorption. Firstly, knowing the sound speed of layers, we calculate, using the Snell-Descartes law (equation 3.1) the angles [θ (l) 0 , θ (l) l-1 ] where θ (l) 0 is the arrival angle detected by the ISM in a homogeneous medium for the layer l. Secondly, we calculate the distance traveled by the wave as follows :

D l = h s cos(θ (l) 0 ) + h r cos(θ (l) 0 ) + l-1 n=1 2h n cos(θ (l) n ) (6.14)
where h n is the thickness of the n th layer.

The density of water is known and the density of the first layer is found in the previous section. The calculation of the remaining buried layers density is done in iterative way. To generalise this method for an arbitrary number of layer l we suppose that we have applied this algorithm to find the densities of the n = l -1 layers. So, in the third step, the impedances of l -1 layers are calculated by the following relation :

Z l-1 = ρ l-1 c l-1 cos(θ (l) l-1 ) (6.15)
The fourth step of the algorithm is the computation of the transmission and reflection coefficients, of a round trip, to each passage in the interface that separates layer n and layer n + 1, from n = 1 to n = l -1 using the following relations :

R n = Z n+1 -Z n Z n+1 + Z n , (6.16)
T n = 1 + R n and T n = -T n . (6.17)

For the l th layer, the reflection coefficient is written as follows :

R l = B l D l A s l-1 n=1 (T n T n ) . (6.18)
Thus, the expression of impedance is obtained as follows :

Z l = Z l-1 1 + R l 1 -R l . (6.19)
Hence, the density is written as :

ρ l = Z l cos(θ l l ) c l . ( 6
.20)

As an example, we illustrate the wave propagation with the reflection and transmission coefficient in each interface to find the density of the fourth layer (l = 4) in figure 6.9.

Case with attenuation

As for the waves propagation where we take into account the absorption in layers, the equation 6.18 of the reflection coefficient is used. This can be done by multiplying the amplitude of the source by the term e -αndn where α n is the absorption coefficient in the layer n and d n is the distance traveled by the wave in the layer n. ESTIMATION

Water (c 0 ) Layer 1 (c 1 , ρ 1 ) Layer 2 (c 2 , ρ 2 ) Layer 3 (c 3 , ρ 3 ) Basement (c 4 , ρ 4 ) x z Source (x s = 0, z s = 0),Amplitude A Receiver z H x H h s h r = h s -z H h 1 h 2 Figure 6
.9 -Propagation of a wave through three layers to find the density of the fourth layer To simplify the calculation, we explain the details step by step, starting from the density of the first layer to generalize the algorithm to an arbitrary number of layers.

To calculate the density of the second layer the equation of direct reflection coefficient R 2 is used as follows :

R 2 = Z 2 -Z 1 Z 2 + Z 1 = B 2K D 2K A s T 1K T 1K exp(θ (2) 1K d 1K ) (6.21)
where Z 1 and Z 2 are respectively the impedance of the first and second layer, B 2K is the amplitude of the reflection waves in the layer 2 recorded on the signal of hydrophone K.

D 2 = h s cos(θ (2) 0K ) + h r cos(θ (2) 0K ) + 2h 1 cos(θ (2) 1K ) , (6.22) d 1K
is the distance taveled by the wave in layer 1 :

d 1K = 2h 1 cos(θ (
= -sin(2φ 1 )(R 2 01 R 12 -R 12 ) (1 + R 12 R 23 cos(2φ 1 )) 2 + (R 12 R 01 sin(2φ 1 )) 2 .
(6.31)

The Z 2 value can be find by solving the non-linear system :

a -a = 0 b -b = 0 (6.32)
To simplify the calculation of the optimization system, we simplify the equations (see appendix B). Then, this system is solved by the fsolve Matlab function with Gauss-Newton algorithm. Hence, the density expression of the second layer is given by :

ρ 2 = Z 2 cos θ (2) 2 c 2 (6.33)
The next step in this algorithm, after finding the density of the second layers, is to estimate the impedance of the following layer in an iterative way by assuming that the first, second and third layers in this algorithm are respectively the layer 1 , the layer 2 and the layer 3 in the stratified medium. The results of this estimation are presented in the section 6.4.

Results

To calculate the density and the absorption coefficient in synthetic data, the two methods (AM and MRM) are used. Then the best one is tested on the tank data.

Synthetic data

To test the efficiency of the proposed algorithms, the acquisition geometry used is the configuration 2 (see subsection 2.2.4) and the seafloor is seabed 2. The geoacoustic parameters are first estimated with MRM algorithm and then with AM algorithm.

(a) The method of multiple reflections

The density of the first layer is calculated using equation (6.13). In applying the MRM algorithm, first, we suppose that we have only three layers which are respectively water, layer 1 and layer 2 then to estimate the density of the third layer we suppose that the three layers are respectively layer 1, layer 2 and layer 3 Results presented in table 6.2 show that the method of multiple reflections gives a correct estimation of the density of layer 2 but not of layer 3 as the error made with the second layer affects the next layer by the fact that the algorithm is very sensitive to small errors. So, the results are not convincing to identify a stratified medium. Indeed, this method requires a careful implementation to reduce its very high risk to digital errors (rounded errors).

RESULTS

139

(b) The amplitude method with and without attenuation

To see the effect of absorption on the amplitude method we perform two tests, initially neglecting absorption and then taking it into account.

Estimated density without taking into account the absorption

Even with a weak absorption, the results of the estimated densities (table 6.3) are affected by this hypothesis. We find out that the theoretical values with the amplitude method are in accordance with the estimated results.

Estimated density taking into account the absorption

The results (table 6.4) show that the densities estimated are very close to the theoretical ones, thus validating the proposed method. But, we note that there is a small shift between the values found and the theoretical ones which is probably due to the uncertainties of input parameters like the arrival angle or to the transmission of the uncertainties from one layer to the others by the fact that we use an iterative algorithm. close to the ground truth in the densities but not for the absorption coefficient.

However, these promising results might be difficult to obtain on data acquired at sea because they are based on an accurate estimation of the amplitudes of the reflected signals. In real experiments, many factors might affect these amplitudes : source directivity, inhomogeneities, presence of very thin layers,. . . Consequently, the amplitudes do not take into account only the impedance contrast but other phenomenons that will degrade the density estimation. The proposed extension of ISM can nevertheless be applied and the densities found could be used as initial starting values for more robust inversion algorithms.

Conclusions

The main objective of this thesis is the estimation of the seabed geoacoustic parameters by improvement and extension of the image source method (ISM). Firstly, this is done by improving the detection and estimation algorithm of the original ISM method, eliminating the false alarms and reducing the computation time. Secondly, the method is extended to the detection of other geoacoustic parameters such as shear wave sound speed and the density.

The main contributions of this thesis can be summarized as follows :

-Study of the ISM algorithmically and theoretically by introducing the idea of virtual image sources in order to understand its advantages and its limitations.

-Improvement of the ISM by using Teager-Kaiser energy operator for detecting the travel times and triangulation for the arrival angles instead of the computation of the migration and the semblance maps.

-Improvement of the new time delay detection algorithm, based on Teager-Kaiser energy operator, by increasing the time resolution using an interpolation and a filtering using the Savitzky-Golay filter followed by the correction of the detection matrix.

-Combination of the Teager-Kaiser energy operator with time-frequency representation to increase the efficiency of the method in very noisy signals environments.

-Development of a new algorithm based on travel time to estimate the true shape of the antenna before inversion.

-Use of the sound speed profile and the arrival angle, previously calculated by the ISM, to estimate the sound speed of shear wave, the densities and the absorption coefficients in sedimentary layers.

These contributions are detailed in the following.

Conclusions

Study of the ISM algorithmically and theoretically by introducing the idea of virtual image sources in order to understand its advantages and its limitations.

To meet the above objectives, a characterization of the seabed geoacoustic parameters is based on modeling the reflection of the wave emitted by a source on a layered medium. Each partial reflection of the wave on the interfaces of the layered medium is represented by an image source whose spatial position is directly related to the speed of the acoustic waves in the layers and their thicknesses.

Both the structure of the seabed sediments which allows us to assume a stratified, homogeneous and parallel fluid medium, and the bistatic configuration which avoids the critical angles, enables us to support the use of such modeling. The information needed to locate the image sources is the source-receiver travel time and the arrival angles on the receiver. Through this modeling, an inversion of the parameters of thickness and sound speed of the layered system originating from the spatial positions of the image sources is made recursively by using the Snell-Descartes law (Chapter 3.3.2).

By locating the image sources in a homogeneous-assumed medium, an error is added to the model. It is due to the choice of the equivalent of the virtual image sources detected in a homogeneous medium and which do not have the same positions in space for each hydrophone. In addition, since the process of determining the sound speed and thicknesses is recursive, it is essential to see whether an error on a layer is propagated and amplified during the course of the algorithm or it is attenuated and how this is done. These errors lead to uncertainty in the inversion result. To make the method applicable, a study has been made, which allows us to retain this.

The choice of the equivalent of the virtual image sources and the median of the antenna have a very little impact on the final inversion results whatever the distance traveled by the wave on the configuration selected. Although this error increases with depth, it remains acceptable. The length of antenna selected to perform the inversion plays a role in the uncertainty of results. The more hydrophones, the more considerable the error is. The axial (travel time) or radial (arrival angle) errors on the positions of the equivalent image sources have far greater effects than the error model. However, the error in time or arrival angle on a layer propagates only to the layer below. The influence of the length of array of hydrophones is the inverse of the error model. The effect of the array length on the error of input parameters has an opposite impact on the model error (choice of equivalent image sources). The longer the antenna, the better results we get. Finally, since the error of input parameters is greater than the model error, the antenna length must be high enough to reduce input parameters error but not so high to reduce the effect of the model error. Thus, 15 or 20 hydrophones is a good compromise between the two errors.

Improvement of the ISM by using a Teager-Kaiser energy operator for detecting the travel times and triangulation for the arrival angles instead of the computation of of the migration and the semblance maps. Since the image sources are reconstructed without refraction, the offsets of the virtual image sources with respect to the vertical of the real source contain the information for determining the sound speed and the layer thickness (SSP). Therefore, the problem is to obtain the positions of these images. This problem of localizing the sources is solved in the original method, by array processing. This processing is based on a signal migration process combined with the use of the semblance function. The problems encountered with the original method are the high computational load of the construction of the migration and the semblance maps and the difficulty in finding the exact position of the virtual image sources in the task generated by the model. Assuming that localizing the virtual image sources allows us to find the two input parameters for the inversion, we have searched for these two parameters in a faster and more accurate manner. To do so, a solution based on the Teager-Kaiser energy operator is proposed.

We have exploited the potential of this quadratic operator to detect the travel time of the image sources directly from the pre-filtered recorded signals (Chapter 4.3). This pre-filtering, which improves the time resolution, is obtained either by cross-correlating the emitted signal with the received one or by using a Wiener deconvolution. The peaks of the output of the energy operator correspond to arrival times but can also be attributed to noise. To overcome this problem, a binary decision strategy based on a threshold is used to isolate the most prominent peaks that are considered as the reflection of the waves on the surface. More precisely, this threshold is introduced to minimize the missing of true peaks, while keeping the number of false detected peaks, caused by noise, within a reasonable limit.

Finally, to obtain the arrival angle required for the inversion algorithm, the intersections located near the nadir of the sources of the pairs of distances d ij and d ik between image source i and hydrophones j and k are first computed. This produces N h (N h -1)/2 locations where N h is the number of hydrophones used. The image source position is taken as the median of all these locations. Equivalently, we compute the median of the hydrophone positions. Finally, to apply the inversion algorithm, from these two positions, the travel time t (i) , and the arrival angle θ (i) for image source i are computed.

Improvement of the new time delay detection algorithm, based on Teager-Kaiser energy operator, by increasing the time resolution using an interpolation and a filtering using the Savitzky-Golay filter followed by the correction of the detection matrix. Two improvements are proposed :

Resolution improvement : the signal is interpolated in order to increase the accuracy of the detected travel time. However, this processing has an impact on noisy signals by introducing false peaks that prevent us to have a reliable automatic detection strategy. To address this problem, a smoothing strategy based on the Savitzky-Golay polynomial filter is used. This filter performs a local polynomial regression to determine the smoothed value for each data point. This smoothing is superior to the moving average because it preserves the characteristics of the data such as the height and width of the peak, which are usually washed by moving average.

Correction of the detection matrix : Under the effect of some disturbances, such as seabed noise, shear waves or differences in the spatial positions of hydrophones, etc, some peaks detected by a hydrophone can not be detected by another. The missing peak, in the detection matrix, for a hydrophone shifts all the values of the associated matrix line. In other words, the detection matrix algorithm replaces the missing values in the matrix by the immediate neighboring values, which makes the matrix completely false. Thus, we have studied the effect of this error by removing some peaks detected on a simulated signal detection matrix, with known ground truth. We have compared the inversion results with the theoretical values. This study allowed us to conclude that the inversion results are still not very far from the ground truth and this is attributed to the choice of the median that is robust to outliers. However, a correction algorithm is imperative in order to correct such a defect created by a missing value. To overcome this problem an automatic detection algorithm has been developed. The correction algorithm has been validated by comparing the theoretical detection matrices and the corrected ones.

Combination of the Teager-Kaiser energy operator with time-frequency representation to increase the efficiency of the method in very noisy signals environments.

To improve the travel time detection accuracy, we noticed that, even though interesting inversion results are obtained, the performances of the inversion method can be affected in very noisy environments. Due to the sensitivity of the Teager-Kaiser energy operator to high-level noise, spurious peaks can be produced, which provide wrong time delay estimates. This is also the case in the presence of impulsive noise as encountered in underwater acoustics. As a solution to this problem, the detection of peaks is done in time-scale or in time-frequency domains where the Teager-Kaiser energy operator is applied to the coefficients at each scale or in each line and each column of spectrogram presentation. This is essentially motivated by the fact that sub-band signals derived, for example, from wavelet decomposition are narrow band where the Teager-Kaiser energy operator is more efficient for peaks detection than on the original signal which can be wide-band. A comparison between the theoretical travel time and the one detected by using different approaches for the synthetic data is done with different SNR values and 15 hydrophones. TK-CWT is effective to detect the desired peaks, even for low SNR (-9 dB), and it reduces the number of false alarms. The computational load of the TK-CWT is lower than that of the two other methods. However, the results of the TK-CWT are conditioned by the choice of the adequate pre-determined mother wavelet and the scale parameter. TK-Spectrogram and TK-ST are efficient and perform similarly for SNR above -4 dB but the TK-ST requires slightly more computational time than the TK-Spectrogram. The main advantage of the TK-ST over the TK-CWT and the TK-Spectrogram, is that it does not require any input parameter. Compared to the TK-CWT, the TK-Spectrogram requires only one parameter, the size of the analyzing sliding window. Finally, for the three methods, the threshold for the binary mask extraction is computed automatically from the data.

Developing a new correcting algorithm based in travel time to estimate the true shape of the antenna before inversion.

The ISM is very sensitive to the acquisition geometry, particularly the geometry of the receiving antenna (Chapter 6.1). The linearity of the antenna is often assumed throughout the experimentation and is not necessarily checked in the case of a long, flexible antenna and towed in a seismic flute type and in the case of an experiment in the sea such as the one used to support this study. Because of the sea currents, the horizontal or vertical antenna may be deformed. Similarly, a dip in the structure may imply that the antenna is not perpendicular or parallel to the seabed. We make no assumption about the shape of the antenna. The correction is reduced to an optimization problem. To estimate the true shape of the antenna, the key step is to use the Teager-Kaiser energy operator to detect the travel time of the direct path and the first reflection recorded by the sensor network. This key step is an alternative to two correction methods with the advantage of fast computation and easy implementation. The first method consists in solving the nonlinear system of the triangulation problem. The second one is to optimize a function of the squared errors Use of the sound speed profile and the arrival angle, previously calculated by the ISM, to estimate the sound speed of shear wave, the densities and the absorption coefficients in sedimentary layers.

For a more complete characterization of seabeds we have extended the ISM to estimate the shear wave of sound speed, the densities and the attenuation coefficients in sedimentary layers using the previously calculated parameters by this method namely, the sound speed, thickness of layers and the arrival angle. The shear wave in a single layer is well detected in synthetic data. Thus, in this case an adaptation of the ISM algorithm provides an accurate estimation of the sound speed. The two analytical methods are : the amplitude method and the multiple reflection method have been used to go back to the two geoacoustic parameters (density and attenuation). The first method uses the data issued from all the hydrophones of the reception antenna and it is based on the amplitudes of direct reflections of the waves on the layers. The second method uses only one hydrophone and it is based on the multiple reflections of the waves on the layers. These two methods, based on iterative estimations of geoacoustic parameters, have been tested on simulated data (vertical and horizontal antenna). The results obtained by the first method are more reliable compared to the theoretical values, and less sensitive to errors that propagate from one layer to the following. Nevertheless, the detection of the amplitudes of reflected signals remains a delicate step. The results of the second method are less convincing. In fact ; this latter requires a careful implementation to reduce its very high risk to digital errors (rounded errors). Overall, the obtained results ; especially with the first method show that a complete characterization of seabed or a classification of sedimentary layers is geologically possible. However, tests on real data are needed with a comparison to field truth to confirm the relevance of the proposed technique associated with ISM. That is why the signals of tank experience have been used. The results are close to the ground truth in the densities but not for the absorption coefficient.

Perspectives

The possible prospects in the direct continuation of this thesis are concerned with two important objectives.

To complete the research track already started, in signal processing contributions, applying the TK-TF approaches on actual data in the sea, SCARAB for example, is a natural and logical continuation of the ideas. The arrival time may advantageously be improved by the application of these approaches and especially for detecting deeper layers in the seafloor than before.

To achieve further work in acoustics, the idea would be to enrich the estimation of stratified medium. First, it is useful to complete the extension of the ISM to find the shear waves in the remaining layers. Then, our estimation technique of density can be used for marine data. For this, the source directivity must be taken into account.

In a mid-term perspective, one important task would be to combine the different versions of ISM. Actually, in recent years, the development of the image source method have provided different methods to localize the image sources. These techniques are the method of migration and semblance maps (ISM-MS), the method proposed in this thesis (ISM-TK), and the approach based on Bayesian inversion processes (ISM-Bayes). Each of these approaches have advantages. By producing a map of the seafloor, ISM-MS can be used for geoacoustic inversion but also for evaluation the seafloor roughness or the volume heterogeneities. ISM-Bayes increases the computation time but provides an estimation of the posterior probability density for the estimated parameters which is useful to quantify the inversion quality but also the heterogeneity degree of a seafloor. ISM-TK provides a very fast and reliable inversion scheme. Moreover, it is rather fully automatic. So, it could be a good candidate for real-time application. Combining these methods can provide better and more accurate results in real-time with optimized geometry.

In this thesis our investigation was limited to shallow water (< 200 m). Seafloor characterization in very shallow water area (< 15 m) is a challenging problem and the ISM is a good candidate tool for such investigation.

Finally, today, the classical inverse methods need huge computing time to estimate the geoacoustic parameters accurately as they do not have a priori information. The image source method with the advantage of fast calculation can provide an initial estimation that can reduce the value intervals to be found. This will ensure much faster calculations.

Appendix A : Migration and semblance maps to localize the virtual image sources positions

The orignal algorithm used for the detection and location of image sources is based on array processing and seismic techniques. It is labelled as ISM-MS. It was developed by S.Pinson [START_REF] Pinson | Caractérisation des fonds marins par la méthode des source images[END_REF] and is described in this Annex. The idea is to migrate the recorded signals in water. They focus on the image sources. Since the focus is not perfect, a semblance map is also compute which, theoretically gives one at image sources location and zero elsewhere. A mask is built on the semblance map where its value is above a given threshold and the maxima of the migration map are picked on this mask (Figure 6.10). This gives the locations of the image sources. Finally, the travel time t (i) and arrival angles θ (i) are computed from these locations and the centre of the array which is its middle position.

Detection and localization of images sources

Migration

The goal is to find the parameters of time and arrival angles between the image sources and the hydrophones. Furthermore, the study is now focused on the reconstruction of the image sources by migrating the signals in a constant speed medium whose speed is that of water.

The preliminary preprocessing phase allows the backscattered signal to gain in temporal resolution and separate each reflection on the interfaces. With the signal quality and since t 0 the instant of the emission of the pulse by the source is known, delays are compensated to find the location of the image sources in the space. Thus, if a source is present at one point, the out of phase signals are sommeront consistently in t 0 . Contrary focusing point r⊥(θ). τ c delay is applied. τ c corresponds to the travel time between the center of the r r c antenna and the point of focus r⊥(θ). Noting d the distance between the center of the r r c antenna and a point r of the image source map (Figure 6.12), the treatment to be applied to the signal S H K (t) of the hydrophone K is : λ is the wavelength in water at the center frequency of the spectrum of the transmitted signal. Therefore :

S H f K (d, θ) = D K (θ)
d ppn = d⊥ 1 + λd⊥ a 2
(6.46)

d dpn = d⊥ 1 -λd⊥ a 2 (6.47)
If λd⊥ a 2 > 1 the last net plane distance is situated in infinity. The distances of the first net plane and the last net plane are drawn in yellow on Figures 6.13 and 6.14 It shows that the deeper image sources are shifted from the vertical but the depth of the field is large enough to allow their detection.

Improved migration feature

The 3 improvements taken into account, the map for displaying the image sources is finally calculated by the function :

I M (d, θ) = 1 N 2 h -N h N h n=1 N h q=1
S H f n (d, θ)S H * f q (d, θ) (6.48)

Semblance

To facilitate an automatic detection of the image source, the semblance function is used. It consists in dividing I M (equation 6.48) by the non-coherent summation the focus signals : The semblance function permits the proper isolation of the small areas where the image sources are located. The problem is that, in a given direction where there is an image, the semblance is 1 on all the duration of the echo, which makes it impossible to know accurately its distance from the center of the antenna. Consequently, I semb (r) is used to define the small areas where there are sources. Then, the sources are accurately located in these areas using I M (d, θ) by searching for maxima premises. Given the discrete nature of the positions of the image sources, even if the resolution methods of the used antennae is limited, the position of the maxima is, in theory, at the exact position of the image. Similarly, in the conventional formation of beamforming, if the main lobe of the antenna directivity has a 60 angle of aperture, this never prevents us from knowing accurately the direction in which we obtain a maximum. However, the positions

I semb (d, θ) = 1 N 2 h -N h

1. 1

 1 Morphology of the marine area (from SEPM Strata website: www.sepmstrata.org/page.aspx). . . . . . . . . . . . . . . . . . . . . . . . 1.2 Depositional environments: the number (1,2,3,4), (5,6,7) and (8,9,10) indicate respectively the origin of marine sediments, the different type of transported and where the sediments are deposited in the sea. (from the University of Washington, Earth and Space science department website: earthweb.ess.washinton.edu/ess-301/Lab2.pdf). . . . . . . . . . . . . . . . 1.3 Example of a core made during the measurement campaign in the region of Bouma [Dev09]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The Sheppard ternary diagram [JR07a]. . . . . . . . . . . . . . . . . . . . 1.5 Sound speed depending on the porosity of the bottom a) on continental shelf b) on the abyssal plains [HB82]. . . . . . . . . . . . . . . . . . . . . 1.6 Density of sediments depending on the grain size [HB82]. . . . . . . . . . . 1.7 Reflection and transmission of plane acoustic wave at a plane interface between two fluid media. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Absorption depending on porosity (left) and on the grain size (right) [HB82]. 1.9 Principle of geoacoustic estimation. . . . . . . . . . . . . . . . . . . . . . . .carleton.edu/eet/seafloor/case-study.html). . . . . . . . . . . . . . . . 1.11 Measurement of marine seabed parameters with monostatic configuration: the sound path is a round trip, a) Side scan sonar b) Multibeam sounder c) Echo sounder [Gui99]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Measurement of seabed parameters by bistatic configuration. . . . . . . . 2.1 Reflection of a wave on a stratified medium in which only the first reflections on the interfaces are shown. . . . . . . . . . . . . . . . . . . . . . . . 2.2 Reflection and transmission of a wave in the case of a layered medium. . . 2.3 Reflection coefficient in the case of multiple fluid interfaces . . . . . . . . . 2.4 Parameters for the calculation of the reflection coefficient of a system of layers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Configuration 1: a vertical antenna. . . . . . . . . . . . . . . . . . . . . . . VIII LIST OF FIGURES 2.6 Configuration 2: a horizontal antenna. . . . . . . . . . . . . . . . . . . . . 2.7 SSP of seabed 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Emitted signal: (left) the signal in temporal domain and (right) in frequency domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Synthetic time signal calculated with numerical evaluation of the Sommerfeld integral in the first hydrophone of the antenna. The numbers 1 to 10 correspond to the reflections on the interfaces of the layered medium. . . . 2.10 Signal synthesized for configuration 2 and seabed 2. The signal correspond to the reflections on the interfaces of the layered medium without direct signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Signal synthesized for configuration 2 and seabed 2 solid. The signal correspond to the reflections on the interfaces of the layered medium without direct signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Experimental tank at ISEN, Brest. . . . . . . . . . . . . . . . . . . . . . . 2.13 Side view of the tank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Top view of the tank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Signal sent in the time and frequency domain . . . . . . . . . . . . . . . . 2.16 Recorded signal in the first hydrophone: 1,2 and 3 are the waves reflected at the different interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Configuration of the tank experience with a horizontal antenna. . . . . . . 2.18 Experiment area of SCARAB measurement. . . . . . . . . . . . . . . . . . 2.19 Experimental signal recorded by the first hydrophone in SCARAB measurement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Reflection of a spherical wave on a simple interface represented by an image source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The refection of a spherical acoustic wave by a stratified fluid bottom. . . 3.3 Modeling of the reflection by the ocean bottom by image sources : each reflection on an interface is replaced by an image source. . . . . . . . . . . 3.4 The equivalent medium of the figure 3.2 by image sources model. . . . . . 3.5 Geometry for the "real" and "virtual" image sources in the case of an increasing SSP [GDCD15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 The two first real and virtual image sources of figure 3.5 for a receiving array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 The VIS of the first and second image sources in the configuration 2 and seabed 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Algorithm of ISM. Green: step 1, Red: step 2. . . . . . . . . . . . . . . . . 3.9 Refracted ray in a stratified medium whose layers thickness and sound speed are known. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES IX 3.10 Inversion of SSP with different number of hydrophones in the case of configuration 1 on seabed 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 1

 1 ISM-MS algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 (a) focusing on the image sources with I BW (r) (in dB). (b): zoom on the image sources. The source is (0,150) and the first image is (0,-150), configuration 1 and seabed 1 [Pin11]. . . . . . . . . . . . . . . . . . . . . . 4.3 Focus on image sources with I semb (d, θ). The black spots correspond to the hydrophone coordinates and the yellow curves represent distances of the first net plane and the last net plane. Configuration 1 and seabed 1 [Pin11]. 4.4 Synthetic signal (left) and its TKEO output (right) in free noise case. . . 4.5 Synthetic signal (left) and its TKEO output (right) with SNR= 3 dB. . . 4.6 Synthetic signal (left) and its TKEO output (right) in noisy case after cross-correlation with the emitted signal. . . . . . . . . . . . . . . . . . . .4.7 The tank recorded signal in Blue and the signal after preprocessing in Red (Left). TKEO output of the recorded signal in Blue and TKEO output of the signal after preprocessing in Red (Right). The solid dots are the travel time detected in each signal. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Synthetic signal preprocessing with matched filtering (Red) or with Wiener deconvolution (Blue). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Synthetic signal after cross-correlation with the emitted signal (Blue). TKEO output of the signal (Red). Peaks detection on TKEO output (Green). True arrival time (Black). . . . . . . . . . . . . . . . . . . . . . . 4.10 The figure presents the threshold T in the different case of the values C in red. The intersection of the blue curve f (m) with the straight black line which presents the value of SNR= 3 dB indicates the value of C to be selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 TKEO outputs of the real tank signals of hydrophones 1 (Blue) and 2 (Red) after preprocessing. The horizontal black line is the calculated threshold. The dots in the curves represent the travel time detected. . . . . . . 4.12 TKEO output of real signals (SCARAB) of hydrophones 1 (Blue) and 2 (Red) after preprocessing. The horizontal black line is the calculated threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X LIST OF FIGURES 4.13 TKEO output of an over-sampled impulse synthetic signal. . . . . . . . . 4.14 In red the TKEO output of an over-sampled impulse synthetic signal of figure 4.13 filtered by SG filter. . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Algorithm of ISM-TK for travel time detection . . . . . . . . . . . . . . . 4.16 Relative errors averaging over 15 hydrophones in each layer between the theoretical time and the detection one with TKEO for different SNR values in the case of synthetic signal without interpolation. . . . . . . . . . . . . 4.17 Relative errors averaging over 15 hydrophones in each layer between the theoretical time and the detection one with TKEO for different SNR values in the case of recorded signals interpolated and filtered with SG filter. . . 4.18 TKEO output of Scarab signals interpolated and filtered for different hydrophones: 1 (Blue),2 (Red),10 (Black) and 15 (Green). On each curves, the solid dots mark the arrival time detected by the algorithm. . . . . . . 4.19 Sound speed error if some peaks are removed. . . . . . . . . . . . . . . . . 4.20 Sound speed error in the case the third peak of the signal is withdrawn on consecutive hydrophones 2, 4, 8, 11, 12 and 14. . . . . . . . . . . . . . . .4.21 Alignment of the TKEO output of hydrophones (1, 5, 10 and 15) of synthetic signals on the first peak. . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Alignment of the TKEO output of interpolated and filtered (SG) SCARAB signals. TKEO outputs of signals of hydrophone 1 (Blue), hydrophone 2 (Red) and hydrophone 15 (Green). . . . . . . . . . . . . . . . . . . . . . . 4.23 Algorithm for detection of missing peaks. . . . . . . . . . . . . . . . . . . 4.24 effects on SSP inversion of the various correction strategies (A to C from top to bottom) for the false negative problem computed on simulated signals.100 4.25 Received signal y(t) on the first hydrophone computed with configuration 2 and seabed 2 (see section 2.2.4) with SNR set to -9 dB. . . . . . . . . . 4.26 Spectrogram of the signal shown in figure 4.25 without (left) and with (right) noise after cross-correlation with the emitted signal. . . . . . . . . 4.27 The scalogram of the received signal y(t) shown in figure 4.25 on the first hydrophone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Plot of E(a) of the different sub-band (scales) a extracted from the scalogram presented in figure 4.27. . . . . . . . . . . . . . . . . . . . . . . . . . 4.29 ST of the noisy signal (SNR=-4 dB) after cross-correlation with the emitted signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.30 TK-Spectrogram of the signal y(t) with SNR=-3 dB. . . . . . . . . . . . . 4.31 TKEO applied to the received signal y(t) with SNR = -9 dB after crosscorrelation with the emitted signal (blue) to its component of scale a = 8 (Red). TKEO applied to the received signal y(t) without noise (black). The dots present the travel time detected. . . . . . . . . . . . . . . . . . . 4.32 Extracted binary mask M T (t, f ) . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES XI 4.33 Relative errors (in percentage) of time-delays estimated by TK-Spectrogram as function of SNR and number of points of the window. . . 4.34 Relative errors (in percentage) of travel time estimated by TK-CWT as a function of Scale and SNR. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.35 Relative errors (in percentage) of time-delays estimated by TK-ST as function of SNR and in a different scales. . . . . . . . . . . . . . . . . . . . . . 5.1 The triangulation to find the coordinates of the point M . . . . . . . . . 5.2 (a) Location of the equivalent image sources relative to different hydrophone positions, with median position labeled by the triangle. (b) Blow up of the rectangle on the left figure. The × label the found positions by the algorithm and the square labels the median of these positions. . . . . . 5.3 Representation of the uncertainty in the sound speed in the case of radius error = 1 m with configuration 1 and seabed 1. . . . . . . . . . . . . . . .

  of the shape of the vertical antenna in SCARAB experiment. Black dots are the hydrophone positions found after correction. Black * are the hydrophone positions before correction. . . . . . . . . . . . . . . . 5.8 Signal received after Wiener deconvolution. 1: Direct 2: first reflection 3: reflection between the two plates. 4: reflection between the second plate and the seafloor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Mode conversion inside a sediment made of one elastic layer above a fluid basement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Example of signals received in the case of a sediment made of a fluid (blue) or an elastic (red) layer over a basement. . . . . . . . . . . . . . . . . . . . 6.3 Detection of arrival time directly on the signals with TKEO in the case of solid layer. TKEO output of this signal (Red). Arrival times detected on the signal (Blue dots). Arrival times detected in the case of fluid layer before (Green dots). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Sketch of the location of image sources in a homogeneous medium: SI 1 ,SI 2 and SI 3 M in the case of fluid layer, SI 1 ,SI 2 and SI 3 S in the case of solid layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XII LIST OF FIGURES 6.5 The equivalent model to inverse the sound speed c s of shear wave in the first layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Strategy to find the density with the ISM . . . . . . . . . . . . . . . . . . 6.7 Reflection of a wave on a stratified medium in which only the first reflections on the interfaces are shown. . . . . . . . . . . . . . . . . . . . . . . . 6.8 Reflection coefficient as a function of incident of angle with (Blue) and without (Red) absorption for the case of configuration 2 and seabed 2 (subsection.2.2.4) simulation. . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Propagation of a wave through three layers to find the density of the fourth layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 ISM-MS algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 (a) focusing on the image sources with I BW (r) (in dB). (b): zoom on the image sources. The source is (0,150) and the first image is (0, -150) . . . . 6.12 Explanatory diagram of the variables used. . . . . . . . . . . . . . . . . . 6.13 Focusing on image sources with I M (r) (dB). The black spots correspond to the contact coordinates of the hydrophones, the yellow curves represent the first net plan and last net plan distances and the red rectangle marks the area used for the calculation of the function I BW (r). In (0, 150 m), the true source can be seen and (0, -150 m) the first image, the background water being at z = 0 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 Focus on image sources with I semb (d, θ). The black spots correspond to the hydrophones coordinates and the yellow curves represent distances of the first net plane and the last net plane. . . . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -Morphology of the marine area (from SEPM Strata website : www.sepmstrata.org/page.aspx).

Figure 1

 1 Figure 1.2 -Depositional environments : the number (1,2,3,4), (5,6,7) and (8,9,10) indicate respectively the origin of marine sediments, the different type of transported and where the sediments are deposited in the sea. (from the University of Washington, Earth and Space science department website : earthweb.ess.washinton.edu/ess-301/Lab2.pdf).

Figure 1

 1 Figure 1.3 -Example of a core made during the measurement campaign in the region of Bouma [Dev09].

Figure 1

 1 Figure 1.4 -The Sheppard ternary diagram [JR07a].

Figure 1

 1 Figure 1.5 -Sound speed depending on the porosity of the bottom a) on continental shelf b) on the abyssal plains [HB82].

Figure 1 . 6 -

 16 Figure 1.6 -Density of sediments depending on the grain size [HB82].

Figure 1

 1 Figure 1.8 -Absorption depending on porosity (left) and on the grain size (right) [HB82].

Figure 1

 1 Figure 1.9 -Principle of geoacoustic estimation.

Figure 1 .

 1 Figure 1.10 -Seismic and sediment sounders (from website : serc.carleton.edu/eet/seafloor/casestudy.html).

Figure 1 .

 1 Figure 1.11 -Measurement of marine seabed parameters with monostatic configuration : the sound path is a round trip, a) Side scan sonar b) Multibeam sounder c) Echo sounder [Gui99].

Figure 1 .

 1 Figure 1.12 -Measurement of seabed parameters by bistatic configuration.

  Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 Synthetic signals . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.1 Reflection of an acoustic spherical wave on a layered seafloor . 32 2.2.2 Plane wave reflection coefficient of a complex seafloor . . . . . 36 2.2.3 Other aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.2.4 Geometric configuration used . . . . . . . . . . . . . . . . . . . 40 2.2.5 Geoacoustic configuration used . . . . . . . . . . . . . . . . . . 42 2.2.6 Examples of signals simulated . . . . . . . . . . . . . . . . . . . 45 2.3 Real signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.3.1 Tank experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.3.2 SCARAB 98 experiment signals . . . . . . . . . . . . . . . . . 51 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Figure 2 . 2 -

 22 Figure 2.2 -Reflection and transmission of a wave in the case of a layered medium.

  Figure 2.5 -Configuration 1 : a vertical antenna.

  Figure 2.6 -Configuration 2 : a horizontal antenna.

Figure 2

 2 Figure 2.7 -SSP of seabed 1.

Figure 2

 2 Figure 2.8 -Emitted signal : (left) the signal in temporal domain and (right) in frequency domain.

Figure 2

 2 Figure 2.9 -Synthetic time signal calculated with numerical evaluation of the Sommerfeld integral in the first hydrophone of the antenna. The numbers 1 to 10 correspond to the reflections on the interfaces of the layered medium.

Figure 2 .

 2 Figure 2.10 -Signal synthesized for configuration 2 and seabed 2. The signal correspond to the reflections on the interfaces of the layered medium without direct signal.

Figure 2 .

 2 Figure 2.11 -Signal synthesized for configuration 2 and seabed 2 solid. The signal correspond to the reflections on the interfaces of the layered medium without direct signal.

Figure 2 .

 2 Figure 2.12 -Experimental tank at ISEN, Brest.

Figure 2 .

 2 Figure 2.13 -Side view of the tank.

Figure 2 .

 2 Figure 2.14 -Top view of the tank.

Figure 2 .

 2 Figure 2.15 -Signal sent in the time and frequency domain

Water

  Figure2.17 -Configuration of the tank experience with a horizontal antenna.

Figure 2 .

 2 Figure 2.18 -Experiment area of SCARAB measurement.

Figure 2 .

 2 Figure 2.19 -Experimental signal recorded by the first hydrophone in SCARAB measurement.

Figure 3

 3 Figure 3.1 -Reflection of a spherical wave on a simple interface represented by an image source.

Figure 3

 3 Figure 3.2 -The refection of a spherical acoustic wave by a stratified fluid bottom.

Figure 3

 3 Figure 3.3 -Modeling of the reflection by the ocean bottom by image sources : each reflection on an interface is replaced by an image source.

2h 2 Figure 3

 23 Figure 3.4 -The equivalent medium of the figure 3.2 by image sources model.

Figure 3

 3 Figure 3.6 -The two first real and virtual image sources of figure 3.5 for a receiving array.

Figure 3

 3 Figure 3.7 -The VIS of the first and second image sources in the configuration 2 and seabed 2.

Figure 3 . 8 -

 38 Figure 3.8 -Algorithm of ISM. Green : step 1, Red : step 2.

Figure 3

 3 Figure3.9 -Refracted ray in a stratified medium whose layers thickness and sound speed are known.

3. 4 .

 4 LIMITS AND SENSITIVITY OF THE ISM TO BAD INPUT PARAMETERS DETECTION. 63differences between the two configurations, which interests us, are the distance traveled by the wave to reach the hydrophone and the type of antenna (vertical or horizontal).

  (a) Inversion with 20 hydrophones. (b) Inversion with 15 hydrophones. (c) Inversion with 10 hydrophones. (d) Inversion with 5 hydrophones.

Figure 3 .

 3 Figure 3.10 -Inversion of SSP with different number of hydrophones in the case of configuration 1 on seabed 1

  Figures 3.10 shows that the algorithm is more effective with 5 hydrophones than with 10 or 15. We move from 0.05 m/s and 0.09 m/s error in the first and the second layers with 5 hydrophones to 1.5 m/s and 2 m/s in the case of 15 hydrophones. The more antenna length is, the more virtual the image sources are extended, and the less accurate the equivalent image position becomes. To see more clearly the effect of antenna length, we plot the uncertainties values curve for different layers with different numbers of hydrophones for configuration 1 and seabed 1 (Figures 3.11).

Figure 3 .

 3 Figure 3.11 -Uncertainties of sound speed (m/s) in each layer in the case of 5,10, 15 and 20 hydrophones for configuration 1 and seabed 1.

Figure 3 .

 3 Figure3.12 -The effect of axial error detection on the image source position.

Figure 3 .

 3 Figure 3.13 -The effect of axial error detection on the sound speed profile

3. 4 .Figure 3

 43 Figure 3.14 -The effect of radial error detection on the image source position

Figure 3 .

 3 Figure 3.16 -Error on sound speed in configuration 1 and seabed 1, in the case of axial and radial position errors in the third layer, with 5,10, 15 and 20 hydrophones

Figure 4

 4 Figure 4.1 -ISM-MS algorithm

Figure 4

 4 Figure 4.2 -(a) focusing on the image sources with I BW (r) (in dB). (b) : zoom on the image sources. The source is (0,150) and the first image is (0,-150), configuration 1 and seabed 1 [Pin11].

Figure 4 . 3 -

 43 Figure 4.3 -Focus on image sources with I semb (d, θ). The black spots correspond to the hydrophone coordinates and the yellow curves represent distances of the first net plane and the last net plane. Configuration 1 and seabed 1 [Pin11].

4. 3 .

 3 NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 81

Figure 4 . 4 -

 44 Figure 4.4 -Synthetic signal (left) and its TKEO output (right) in free noise case.

Figure 4 . 5 -

 45 Figure 4.5 -Synthetic signal (left) and its TKEO output (right) with SNR= 3 dB.

Figure 4 .

 4 Figure 4.6 shows the importance of the preprocessing step before the application of the

Figure 4 . 6 -

 46 Figure 4.6 -Synthetic signal (left) and its TKEO output (right) in noisy case after crosscorrelation with the emitted signal.

Figure 4

 4 Figure 4.7 -The tank recorded signal in Blue and the signal after preprocessing in Red (Left). TKEO output of the recorded signal in Blue and TKEO output of the signal after preprocessing in Red (Right). The solid dots are the travel time detected in each signal.

Figure 4 . 8 -

 48 Figure 4.8 -Synthetic signal preprocessing with matched filtering (Red) or with Wiener deconvolution (Blue).

Figure 4

 4 Figure 4.9 -Synthetic signal after cross-correlation with the emitted signal (Blue). TKEO output of the signal (Red). Peaks detection on TKEO output (Green). True arrival time (Black).

Figure 4 .

 4 Figure 4.10 -The figure presents the threshold T in the different case of the values C in red. The intersection of the blue curve f (m) with the straight black line which presents the value of SNR= 3 dB indicates the value of C to be selected.

Figure 4 .

 4 Figure4.11 shows that the threshold value is perfectly selected, by signals averaging, compared to the ground truth for the two hydrophones 1 and 2 which have been chosen for illustration in this analysis.

Figure 4 .

 4 Figure 4.11 -TKEO outputs of the real tank signals of hydrophones 1 (Blue) and 2 (Red) after preprocessing. The horizontal black line is the calculated threshold. The dots in the curves represent the travel time detected.

Figure 4 .

 4 Figure 4.12 -TKEO output of real signals (SCARAB) of hydrophones 1 (Blue) and 2 (Red) after preprocessing. The horizontal black line is the calculated threshold

4. 3 .

 3 NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY OPERATOR 87

Figure 4 .

 4 Figure 4.13 -TKEO output of an over-sampled impulse synthetic signal.

Figure 4 .

 4 Figure 4.14 -In red the TKEO output of an over-sampled impulse synthetic signal of figure 4.13 filtered by SG filter.

Figure 4 .

 4 Figure 4.15 -Algorithm of ISM-TK for travel time detection

Figure 4 .

 4 Figure 4.16 -Relative errors averaging over 15 hydrophones in each layer between the theoretical time and the detection one with TKEO for different SNR values in the case of synthetic signal without interpolation.

Figure 4 .

 4 Figure 4.17 -Relative errors averaging over 15 hydrophones in each layer between the theoretical time and the detection one with TKEO for different SNR values in the case of recorded signals interpolated and filtered with SG filter.

Figure 4 .

 4 Figure 4.18 -TKEO output of Scarab signals interpolated and filtered for different hydrophones : 1 (Blue),2 (Red),10 (Black) and 15 (Green). On each curves, the solid dots mark the arrival time detected by the algorithm.

  Figure 4.19 -Sound speed error if some peaks are removed.

Figure 4 .

 4 Figure 4.20 -Sound speed error in the case the third peak of the signal is withdrawn on consecutive hydrophones 2, 4, 8, 11, 12 and 14.

Figure 4 .

 4 Figure 4.21 -Alignment of the TKEO output of hydrophones (1, 5, 10 and 15) of synthetic signals on the first peak.

  Figure 4.22 -Alignment of the TKEO output of interpolated and filtered (SG) SCARAB signals. TKEO outputs of signals of hydrophone 1 (Blue), hydrophone 2 (Red) and hydrophone 15 (Green).

Figure 4 .

 4 Figure 4.23 -Algorithm for detection of missing peaks.

Figure 4 .

 4 Figure 4.24 -effects on SSP inversion of the various correction strategies (A to C from top to bottom) for the false negative problem computed on simulated signals.

Figure 4 .

 4 Figure 4.25 -Received signal y(t) on the first hydrophone computed with configuration 2 and seabed 2 (see section 2.2.4) with SNR set to -9 dB.

Figure 4 .

 4 Figure 4.26 -Spectrogram of the signal shown in figure 4.25 without (left) and with (right) noise after cross-correlation with the emitted signal.

  ).

Figure 4 .

 4 Figure 4.27 -The scalogram of the received signal y(t) shown in figure 4.25 on the first hydrophone.

Figure 4 .

 4 Figure 4.28 -Plot of E(a) of the different sub-band (scales) a extracted from the scalogram presented in figure 4.27.

S

  T {data} = S T {signal} + S T {noise}. (4.37) An example of ST applied in the signal y(t) with SNR=-4 dB is shown in figure 4.29.

Figure 4 .

 4 Figure 4.29 -ST of the noisy signal (SNR=-4 dB) after cross-correlation with the emitted signal.

Figure 4 .

 4 Figure 4.30 shows the spectrogram of y(t) with SNR=-3 dB and the associated extracted mask M p (t, f ). The peaks are well evidenced on the mask.

( a )

 a Spectrogram of y(t).(b) Extracted binary mask M p (t, f ).

Figure 4 .

 4 Figure 4.30 -TK-Spectrogram of the signal y(t) with SNR=-3 dB.

Figure 4 .

 4 Figure 4.31 -TKEO applied to the received signal y(t) with SNR = -9 dB after cross-correlation with the emitted signal (blue) to its component of scale a = 8 (Red). TKEO applied to the received signal y(t) without noise (black). The dots present the travel time detected.

Figure 4 .

 4 Figure 4.34 -Relative errors (in percentage) of travel time estimated by TK-CWT as a function of Scale and SNR.

Figure 4 .

 4 Figure 4.35 -Relative errors (in percentage) of time-delays estimated by TK-ST as function of SNR and in a different scales.

5. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Image sources localization for input parameters estimation . 5.2.1 Triangulation principle . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Analysis of the sensitivity . . . . . . . . . . . . . . . . . . . . . 5.2.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Deformation and correction of the array . . . . . . . . . . . . 5.3.1 Triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Optimisation error approach . . . . . . . . . . . . . . . . . . . 5.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Inversion results . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.2 Tank experiments . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.3 Inverse results of real data Scarab 98 . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5

 5 Figure 5.2 -(a) Location of the equivalent image sources relative to different hydrophone positions, with median position labeled by the triangle. (b) Blow up of the rectangle on the left figure. The × label the found positions by the algorithm and the square labels the median of these positions.

Figure 5

 5 Figure 5.3 -Representation of the uncertainty in the sound speed in the case of radius error = 1 m with configuration 1 and seabed 1.

Figure 5

 5 Figure 5.4 -Representation of the distances between the source and the hydrophones and the distances between the first image source and the hydrophones

Figure 5

 5 Figure 5.7 -Correction of the shape of the vertical antenna in SCARAB experiment. Black dots are the hydrophone positions found after correction. Black * are the hydrophone positions before correction.

Figure 5

 5 Figure 5.8 -Signal received after Wiener deconvolution. 1 : Direct 2 : first reflection 3 : reflection between the two plates. 4 : reflection between the second plate and the seafloor.

6. 2 .

 2 SHEAR WAVE SOUND SPEED INVERSION 127 Water (c 0 , ρ 0 ) Layer 1 (c 1 , ρ 1 ) Basement (c 2 , ρ 2 ) P P

Figure 6

 6 Figure 6.1 -Mode conversion inside a sediment made of one elastic layer above a fluid basement.

Figure 6

 6 Figure 6.2 -Example of signals received in the case of a sediment made of a fluid (blue) or an elastic (red) layer over a basement.

Figure 6 . 3 -

 63 Figure 6.3 -Detection of arrival time directly on the signals with TKEO in the case of solid layer. TKEO output of this signal (Red). Arrival times detected on the signal (Blue dots). Arrival times detected in the case of fluid layer before (Green dots).

WaterFigure 6

 6 Figure 6.4 -Sketch of the location of image sources in a homogeneous medium : SI 1 ,SI 2 and SI 3 M in the case of fluid layer, SI 1 ,SI 2 and SI 3 S in the case of solid layer.

Figure 6 . 5 -

 65 Figure 6.5 -The equivalent model to inverse the sound speed c s of shear wave in the first layer.

  -Descartes laws of refraction, the shear sound speed is estimated by :

Figure 6

 6 Figure 6.7 -Reflection of a wave on a stratified medium in which only the first reflections on the interfaces are shown.

Figure 6 . 8 -

 68 Figure 6.8 -Reflection coefficient as a function of incident of angle with (Blue) and without (Red) absorption for the case of configuration 2 and seabed 2 (subsection.2.2.4) simulation.

Figure 6 .

 6 Figure 6.10 -ISM-MS algorithm

  (θ) -τ c (θ)) (6.41) with τ K (θ) = |r r K -r⊥(θ)| c 0 and τ c (θ) = |r r c -r⊥(θ)| c 0 . The multiplication by D K (θ) =|r r K -r⊥(θ)| is the compensation of the spherical geometric divergence. By focusing on the vertical antenna, the sources that are too far apart will not be detected. The image sources are generally close to the vertical but it should be ensured that the depth of the focalization field is large enough to detect the image sources slightly offset from the vertical. The field depth is the area between the distance of the first net plane (d ppn ) and the area of last net plane (d dpn )[START_REF] Pérez | Optique géométrique et ondulatoire avec 200 exercices et problèmes résolus[END_REF] :d ppn = d⊥ -res ppn d⊥ a (6.42)d dpn = d⊥ + res dpn d⊥ a (6.43)where d⊥ is the focusing distance |r r c -r⊥(θ)| and |r r N -r r 1 |sin θ is the apparent scale of the antenna. res ppn and res dpn are respectively the limits of transverse resolutions for the first net plan and the last net plan in the axis of the antenna :

  f n (d, θ)S H * f n (d, θ) (6.51) With this function, the result is the value 1 if all the focused signals are identical in phase and amplitude, and a value close to 0 if they are different from each other (Figure 6.14).

Figure 6 .

 6 Figure 6.14 -Focus on image sources with I semb (d, θ). The black spots correspond to the hydrophones coordinates and the yellow curves represent distances of the first net plane and the last net plane.

  Contents 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Marine sediments . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 1.2.2 Structure of marine sediments . . . . . . . . . . . . . . . . . . . 8 1.3 Geological property of marine sediments . . . . . . . . . . . . 10 1.3.1 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 Porosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3.3 Gradients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 Geoacoustic parameters . . . . . . . . . . . . . . . . . . . . . . 14 1.4.1 Sound speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.2 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.3 Acoustic impedance . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.4 Attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.5 Characterization of seabed by acoustic methods . . . . . . . . 21 1.5.1 Why are the acoustic methods used ? . . . . . . . . . . . . . . 21 1.5.2 Geoacoustic systems . . . . . . . . . . . . . . . . . . . . . . . . 21 1.5.3 Geoacoustic inversion methods . . . . . . . . . . . . . . . . . . 25 1.5.4 The image source method . . . . . . . . . . . . . . . . . . . . . 28 1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Table 1 .

 1

1 -Wentworth classification of sediment

[START_REF] Wentworth | A scale of grade and class terms for clastic sediments[END_REF]

.

  , LMR + 92]. Table 1.2 shows some minimum and maximum values of sound speed and densities in different types of seabed sediments.

	Sediment Type Sound speed (m/s) Density (kg/m 3 )
		Min	Max	Min	Max
	Coarse sand	1780	1880	2030	2080
	Medium sand	1830	1800	2000	2040
	Fine sand	1645	1700	1920	1970
	Very fine sand	1680	1710	1880	1940
	Silty sand	1540	1660	1680	1800
	Silty sand	1490	1645	1660	1790
	Silt	1495	1620	1650	1740
	Sand-silt-clay	1480	1590	1450	1620
	Clay loam	1450	1550	1430	1480
	Table 1.2 -Minimum and maximum values for sound speed and density measured in sedi-
	ments [Rak09].

Table 1 .

 1 3 -Conversion of the most used units for absorption coefficient where f is the frequency and the c is the sound speed.

	1.3.

  Table 2.1 -Coordinate of hydrophones and their spaced in the plan (x,z).

	Hydrophones	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
	x (m)	200	200	200	200	200	200	200 200 200	200	200	200	200	200	200
	z (m)	11.5 15.5 19.5 23.5 28.5 32.5	39	45	51	55.5 58.5 61.5 65.5 69.5 73.5
	Hyd-spaced(m)	-	4	4	4	5	4	6.5	6	6	4.5	3	3	4	4	4

  .7.

	Layers	Sound speed (m/s) Density Thickness (m)
	Water	1500	1	.
	layer 1	1520	1.1	3
	layer 2	1540	1.2	3.5
	layer 3	1600	1.5	4
	layer 4	1630	1.7	2
	layer 5	1700	1.9	6
	layer 6	1720	2	2
	layer 7	1800	2.5	5
	layer 8	1900	3	3
	layer 9	1920	3.1	2
	Basement	2000	3.6	.

Table

2

.2 -The geoacoustic parameters of the layered medium labeled seabed 1. The absorption coefficient of layers are assumed zero.

Table 2

 2 Thickness Density sound speed c p sound speed c sThe seafloor is composed of 9 fluid layers on a basement, the geoacoustic parameters of these layers are presented in the table 2.4. The absorption is assumed zero in the medium.

	Absorption

.4 -The geoacoustic parameters of the layered medium labeled seabed 3. The absorption coefficient of layers are assumed zero.

Table 3

 3 Table 3.2 -Error of sound speed and thickness in the case of the three different seabed with the geometric configuration 2.

	Layers	Error of Sound speed 10 -2 (%)	Error of thickness 10 -2 (%)
		Seabed 1 Seabed 3 Seabed 4 Seabed 1 Seabed 3 Seabed 4
	layer 1	0.05	0.18	-	2.39	4.5	-
	layer 2	0.08	0.03	0.1	2.41	0.04	3.13
	layer 3	0.34	0.92	0.28	3.94	6.75	3.56
	layer 4	0.29	0.05	-	3.17	2.16	-
	layer 5	0.42	0.42	0.45	3.27	3.2	3.5
	layer 6	0.11	0.24	-	1.92	2.23	-
	layer 7	0.11	0.17	0.06	1.6	0.53	1

.1 -Error of sound speed and thickness in the case of the three different seabed with the geometric configuration 1 by 15 hydrophones with 15 m length.

  OF THE ISM TO BAD INPUT PARAMETERS DETECTION. 65 Layers Error of Sound speed (%) 10 -2 Error of thickness (%) 10 -2 Seabed 1 Seabed 2 Seabed 3 Seabed 1 Seabed 2 Seabed 3 Table 3.3 -Error of sound speed and thickness in the case of the three different seabed with the configuration 1 by 15 hydrophones with 64 m length.

	layer 1	2.6	4.81	-	3.30	8.12	-
	layer 2	4.54	2.15	4.65	5.71	3.31	7.68
	layer 3	10.65	13.8	10.42	17.54	23.73	17.12
	layer 4	12.88	9.89	-	20	15.78	-
	layer 5	17.64	17.7	17.88	28.3	29	29.48
	layer 6	18	23.3	-	25	38.46	-
	layer 7	22.7	15.17	22.62	36	21.52	36.38

  Contents 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Original ISM Algorithm for detection and localization . . . . 76 4.3 New detection algorithm based on the Teager-Kaiser energy operator .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3.1 Time delay estimation problem . . . . . . . . . . . . . . . . . . 79 4.3.2 Definition of the signals energy . . . . . . . . . . . . . . . . . . 80 4.3.3 Teager-Kaiser Energy Operator . . . . . . . . . . . . . . . . . . 80 4.3.4 TKEO frequency Analysis . . . . . . . . . . . . . . . . . . . . . 82 4.3.5 TKEO Application . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.6 The choice of the preprocessing . . . . . . . . . . . . . . . . . 84 4.3.7 Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.8 Resolution improvement . . . . . . . . . . . . . . . . . . . . . 90 4.3.9 Improvement of matrix detection . . . . . . . . . . . . . . . . . 94 4.4 TKEO associate with time-frequency representations . . . . . 101 4.4.1 Time-frequency representations . . . . . . . . . . . . . . . . . . 101 4.4.2 TKEO in time-frequency domain . . . . . . . . . . . . . . . . . 104 4.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

  a simple model of mechanical acoustical system, where the object oscillates, thus creating pressure waves in the surrounding medium. The solution of equation (4.8) consists of the oscillating signal x(t) = A cos(2πf t + φ 0 ), where A is the amplitude, f = 1 2π

	4.3. NEW DETECTION ALGORITHM BASED ON THE TEAGER-KAISER ENERGY
	OPERATOR	79
	Equation 4.8 is
	.8)

  and f (t) do not vary too fast or too much with respect to the average value of f (t)[START_REF] Maragos | On amplitude and frequency demodulation using energy operators[END_REF]. Equation (4.17) shows that the TKEO depends upon the instantaneous amplitude and IF which are attributed to the physics of the signal generation. Given the simplicity of the TKEO and the broad applicability of the AM-FM model in signal processing and communication systems, this operator has found many applications in various engineering domains and particularly via the AM-FM demodulation scheme introduced by Maragos et al., the energy separation algorithm (ESA)[START_REF] Maragos | Energy separation in signal modulations with applications to speech analysis[END_REF].

		t) dt 2	(4.15)
	When applied to an AM-FM signal
	t
	x(t) = a(t) cos 2π	f (τ )dτ	(4.16)
	0
	with time varying amplitude envelope a(t) and instantaneous frequency (IF) f (t), the Ψ
	operator yields
	Ψ[x(t)] ≈ a 2 (t)f 2 (t)	(4.17)
	assuming that a(t)

Table 4

 4

.1 -Matrix detection of travel time in the first case.

Table 5 .

 5 2 -Comparison of the geoacoustic inversion results in the simulated environment before and after the correction of the antenna shape.

	Layer	Type	Theoretical values without correction with correction
	Layer 1	Thickness (m)	2	1.49	2
		Sound speed(m/s)	1650	1482	1650.1
	Layer 2	Thickness (m)	2	1.35	2
		Sound speed (m/s)	1750	1664.1	1750.1

  Table5.4 -Sound speed and thickness values of synthetic data estimated by ISM-TK and ISM-MS methods in the case of configuration 2 and seabed 2.

			values ISM-MS ISM-TK (2 Hy) ISM-TK (15 Hy)
	Layer1	Thickness (m)	2	2.3	2	1.99
		Sound speed (m/s)	1650	1658.7	1650	1650.1
	Layer2	Thickness (m)	4	4	4	3.98
		Sound speed (m/s)	1750	1744.2	1750	1750.1

Table 5 .

 5 8 -Sound speed and thickness values of SCARAB data estimated by ISM-TK, ISM-MS and Holland & Osler methods. ". . . " indicated that the layer has not been identified.

	Holland & Osler		ISM-MS			ISM-TK
			Sound			Sound			Sound
		Depth speed		Depth speed		Depth speed
	Layer	(m)	(m/s) Layer	(m)	(m/s) Layer	(m)	(m/s)
	1	0.5	1502	. . .
	2	1.1	1551	a	1.1	1464	I	1.4	1560
	3	3.3	1516	b	2.8	1496	II	2.99	1517.9
	4	4.8	1527	c	4.3	1606	III	6.4432 1526.8
	5	5.6	1591	d	5.1	1516	. . .
	6	15.1	1555	e	7	1530	IV	14.191 1542.3
				f	14.7	1584

  R 01 R 12 cos(2φ 1 )) + R 12 (cos(2φ 1 ) + R 12 R 01 )) (1 + R 12 R 01 cos(2φ 1 )) 2 + (R 12 R 01 sin(2φ 1 )) 2

		CHAPTER 6. EXTENSION OF ISM TO OTHER GEOACOUSTIC PARAMETERS
				ESTIMATION
	with	a =	(R 01 (1 + ,	(6.30)
	and
			b
			1K ) 2)	(6.23)
	Using equation 6.21 we get :

  Table 6.2 -Estimated results of the method of multiple reflections

		Theoretical density Estimated density
	layer 1	1.8	1.8
	layer 2	2	1.92
	layer 3	2.2	3.5

on simulated data sets used as examples in Chapter 2.

In the next chapters we present the methods used to obtain the travel time and angles directly from the signal, without going through the maps as it was in the original method in order to save time in calculating them and eliminate more the false alarms. and each scale, a, the relative error between the true time delay and the estimated one by the method TK-CWT is shown in figure 4.34. This figure shows that the relative errors

CHAPTER 5. IMAGE SOURCES LOCALIZATION AND INVERSION RESULTS

the original method ISM-MS. In the next chapter, we will see the extension of ISM method to find other parameters to better characterize the seabed.

The theoretical sound speed c s Estimated sound speed c s 700 m/s 707 m/s -first by the fact that we use the estimated parameters c1 and h1 in the calculation instead of the real one. The sensitivity study of this type of error is not yet done. -second by the fact that we use the equivalent image source to find the input parameters which means model error.

In summary, the image source due to shear wave is clearly detected and, with this appropriate model and algorithm presented in this section, shear sound speed can be estimated from the data. These numerical results, which must be confirmed by a deeper study, open the door for a way to invert both compressional and shear properties of sediments with ISM.

Estimation of the density with the ISM

Strategy

The ISM provides two geoacoustic parameters (sound speed and thickness) for one layer from the knowledge of the two input parameters estimated from the data (travel time and arrival angle). Estimating the density requires an additional input. For this, we use the signals amplitude. The method uses the principle of wave propagation in a stratified medium under the Born approximation (subsection 2.2.1) by neglecting the multiple reflections. The knowledge of the sound speed, the thickness and the arrival angles, allows us to estimate the density by measuring the values of the amplitudes of the emitted and reflected signals to compute the reflection coefficient. We are interested by the local reflection coefficient R ij at each interface and not the global reflection coefficient of the stratified medium. The strategy to find the density is presented in figure 6.6. 0K : the arrival angle of the first reflection in the receiver K. θ (l) 0K : the arrival angle of the l th reflection in the receiver K.

All equations are written with 15 hydrophones by replacing each time the recorded amplitudes and the arrival angle. Finally, the optimization of this system enables us to determine the absorption of the first layer α 1 and the density of the second layer ρ 2 . To do so, we use the MATLAB function "fsolve" which enables us to calculate the solution of non linear system.

For the l th layer, with K hydrophones the equation is written as :

The results of this estimation are presented in the section 6.4.

Multiple reflection method (MRM)

As we have seen in the section 2.2.2, to find the global reflection coefficient in three fluid layers we have used two different methods. Thus, two different equations have been find :

-first the equation 2.24, the global reflection coefficient as a function of the local reflection coefficients, the absorption coefficients, and the thicknesses. -second the equation 2.33, the global reflection coefficient as a function of the input impedance, the absorption coefficients, and the thicknesses. In this section, the density is calculated by using the equality of these two equations.

As the density of water is known and the density of the first layer is found in the subsection 6.3.3, the aim is to determine the impedance of the second layer by solving the equality equations. Then we find the impedances of the remaining layer in an iterative way.

We start by equalizing the two equations to find the impedance Z 2 . The calculation is done by separating the real part from the imaginary part. Equation 2.33 is then written :

(6.28)

The equation 2.24 becomes : R = a + ib (6.29) CHAPTER 6. EXTENSION OF ISM TO OTHER GEOACOUSTIC PARAMETERS ESTIMATION

Results of real tank experience

The ISM algorithm is applied on the signals recorded in the sensor after preprocessing and provides the thickness and sound-speed for the two plates. Then, the proposed algorithm leads to an estimation of the density of these two layers.

The results are displayed in Table 6.5. For the first layer made of Teflon, the result is very accurate. For the resin layer, the estimated density is less accurate but still acceptable. The algorithm also provides an estimations of the absorption coefficient but the results are far from classical values and, moreover, can not be compared to any ground truth. Table 6.5 -Geoacoustic parameters of ground truth

Method

Conclusion

In this chapter, we have extended the ISM to estimating the shear wave sound speed, the densities, and the absorption coefficients of sedimentary layers using the parameters estimated by the ISM in the previous chapters, namely the pressure wave sound speeds, the thickness of layers and the arrival angles.

In the first part of this chapter the shear wave in a single layer is well detected in synthetic data. An adaptation to this case of ISM algorithm provides an accurate estimation of the sound speed. This result opens the door to an interesting application of ISM. Indeed, the problem of modeling sediments as fluid or visco-elastic media is still in debate. Application of ISM to data acquired at sea in various sedimentary environments could bring some insights to this debate.

In the second part, two analytical methods are used to estimate the densities and the absorption coefficient : the amplitude method and the multiple reflection method. The obtained results by the AM method are more reliable compared to the theoretical values, and less sensitive to errors that pass on from one layer to another. Nevertheless, practically the detection of the amplitudes of reflected signals remains a delicate step. Overall, the obtained results, especially with the AM, show that a complete characterization of seabed or a classification of sedimentary layers is geologically possible. The method is also tested on real data. The tank experiments reveal the quality of the proposed algorithm with a comparison of true values to the estimated ones. The results are Appendix A : Migration and semblance map to localize the image sources method to what was mentioned in the literature, Kirchhoff migration [START_REF] Claerbout | Downward continuation of moveout corrected seismograms[END_REF], the signals are migrated till the image sources and not till the geological interfaces. This migration of recorded signals allows to focus the source position as well as the image sources corresponding to the reflections on the interfaces in a medium without a sedimentary structure.

As it has shown, the signals from the image sources are only reflections of the signal sent from the source on the interfaces which explains cohesion. Accordingly, the signal received by the hydrophone K and the coordinates r sm is only the sum of pulses with different amplitude factors m sent by all the sources as follows :

where * is the convolution product, g(r s , r r K , t) is the Green function in the temporal domain, η K (t) is an additive noise assumed spatially white and f (t) is the transmitted signal. r s and r r K are respectively the position of the source and the receiver. In this model, the function Green g(r s m , r r K , t) is defined by :

where δ(t) is the Dirac function, D(r s m , r r K ) is the term of geometric divergence and τ (r s m , r r K )) is the travel time between the source and the sensor. In the frequency domain, the expression is written :

with G(r s , r r K , ω) and F (ω) the Fourier transform of g(r s , r r K , t) and f (t). Recall that the coordinates of r s m images are indexed from 0 and the coordinate of the true source has no index. In order to retro propagate the signals at r, the S K spectrum (ω) is multiplied to the inverse of the Green function in a homogeneous medium (c 0 speed)

This operation is performed in the frequency domain because it allows us to obtain very accurate phase shifts in contrast to the application of a delay in the time domain. Indeed, in the time domain, the accuracy of delay is limited by the sampling frequency of the signal.

The retro propagated signals are then filtered by a Gaussian window w(t) which has approximately the same duration as the transmitted signal f (t), then the average of the N h hydrophones is calculated. The result (homogeneous to pressure) is mapped by :

Equivalently, I BW (r) can be written as the sum of the elements of the covariance matrix :

S bn (r, t)S * bq (r, t)dt (6.39)

where S * bq (r, t) is the conjugate of S bq (r, t). In this presentation (Figure 6.11) it appears that there are many artifacts whose amplitudes are sometimes larger than the image sources which are visible as well as the source. The location of these image sources are just the intersections of the circles with the radii |r s 0 -r r n | centered on the n sensor. The values of the amplitudes are summed in these intersection areas to form a spot. The problem with this map is that the artifacts which present isolated circles issued from the real source , which have , in certain cases , an amplitude superior to the image sources. As a result, an automatic detection will be difficult to achieve.

Improving the image

The purpose of this step is to avoid windowing (for removing the oscillation from the pulses) which deteriorates the resolution in distance. Windowing is replaced by demodulating the pulses. For this purpose, the Hilbert transform is used [START_REF] Poularikas | The transforms and applications handbook[END_REF], working with analytical signals. The signal is demodulated in amplitude. The analytical signal S H K (t) of the S K (t) signal is defined by : Appendix A : Migration and semblance map to localize the image sources method

where H(S K (t)) is the Hilbert transform of S K (t).

The second approach achieved to improve the results is to consider the diagonal can improve the antenna directivity function.

Improving the calculation time

The geometric variables used in this section are explained in the diagram. Calculating the I BW (r) function involves a calculation of the phase shift of the recorded signals S K (t) different for all r points in space and thus a significant calculation cost. To decrease this cost, we use the fact that the image sources are supposed to be near the vertical of the actual source. For this, we apply the delay laws of the signals so that the antenna is directed at an angle θ, while focusing on a point of the vertical axis of the actual source r⊥(θ). Thus, for each incidence guidance θ of the antenna, a single phase shift is to be calculated for all r points located on the right incidence θ passing through the center of the r r c antenna (and through the focusing point r⊥(θ)). To focus the signals on the vertical of the source, a lead τ K corresponding to the travel time of the wave between the hydrophone K and the focusing point r⊥(θ) is applied to them. Thus, the phase-shifted signal no longer has the arrival time information for determining the distance between the image source and the antenna. To apply , for this signal, that travel time information from the center of the antenna r r c until the where N 2 h -N h is the number of elements of the double sum and the root is used to obtain a homogeneous pressure result. The result I M (d, θ) result is calculated as a function of d and θ. It is therefore drawn along Straight implications θ through the center of the antenna. To convert it to function r = (x, z), just do the following change of variables :

The result of the I M function (equation 6.48) is shown in Figure 6.13. In this figure, the parasitic circles drawn by the echoes are lower in amplitude about 7 dB to those traced by I BW (equation 6.38) as the elements S H f n (d, θ)S H * f nq (d, θ) with n = q are not taken into account in the double sum of equation 6.48. The attenuation of the circles should normally be much larger because the simulated data are used without adding any noise. However, the numerical calculation of the Sommerfeld integral generates a slight background noise that allows the appearance of these circles when calculating I M . Although the circles are still very attenuated their presence involves the same problems with I BW , that is, it is difficult to make an automatic detection of the image sources from I M . We can see that the deepest image source are shifted vertically (about ten meters for the deepest) but the depth of field is large enough to allow their detection (yellow curves in figure 6.13).

Figure 6.13 -Focusing on image sources with I M (r) (dB). The black spots correspond to the contact coordinates of the hydrophones, the yellow curves represent the first net plan and last net plan distances and the red rectangle marks the area used for the calculation of the function I BW (r). In (0, 150 m), the true source can be seen and (0, -150 m) the first image, the background water being at z = 0 m.

Appendix A : Migration and semblance map to localize the image sources method of the maxima can be distorted by the presence of noise. The resolution limit in the proposed method will intervene mainly in the ability to distinguish very similar images corresponding to thin strata.

Appendix B : Bayesian method to localize the virtual image sources positions

Another way, rather different of the previous one, to obtain the position of the virtual image sources is to use Bayesian inversion [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. The starting point here is to consider that we are seeking the median position (x, z) of virtual image sources of a particular layer i and that t K with K ∈ {1, N h } are the travel times between the hydrophones and this median position. As previously, these travel times are directly converted into distances d K . Note again that the exponents (i) are not displayed for clarity. As the input quantities of the SSP inversion algorithm are the arrival times and the incident angles, the position of the virtual image source is seek in polar coordinates with a new coordinate system centred on the median position of the array. The position (x H (K), z H (K)) of the hydrophones are converted in to polar coordinates (ρ H (K), φ H (K)). The model parameters are then the coordinates of the virtual image source : m = (ρ v , φ v ). The data are the arrival travel times : d = t K with K ∈ {1, N h }. We assume that the hydrophones positions and the water sound-speed are perfectly known. But, within the context of Bayesian inversion, this hypothesis could be relaxed and replaced by the knowledge of prior information. The forward model d = g(m) is simply solved by the following system of equations :

The objective is to obtain the Posterior Probability Density (PPD) expressed in the Bayes theorem :

Once the experiments have been made and the data obtained, P (d) is a constant factor and the probability density function (PDF) P (d/m) can be identified to the likelihood function. If we assume that the uncertainties on the data are Gaussian, independent, and with a common standard deviation σ, this likelihood function can be written as :

where Γ is a normalization factor, the arrival times t cal K are obtained with the forward Appendix B : Bayesian method to localize virtual image source position model, and the arrival times t obs K are the ones obtained with the TK detection algorithm.

At this stage, we do not set any priority to the location of the virtual image source, so, the prior PDF P (m) is set to a constant. We then have set all the necessary pdf. Since the model parameter is only in 2 dimensions, the most easiest way is to define a grid for ρ v and φ v and to compute the PPD at each point of this grid and this for the N s detected image sources. The position of the virtual image source is taken as the maximum of the corresponding PPD. Once again, we therefore have all the input data for the SSP inversion algorithm. Compared to the two other algorithms, this procedure has the advantage of being much more straightforward and of avoiding the computation of maps and the solution of a nonlinear system. Moreover, it gives not only the value of the input data but also a quantification of the uncertainty on these data Appendix C : Simplified equations for the multi-reflection method

The reflection coefficient can be written in two ways :

For any layer P , with p ∈ [1, n -1], we have :

To simplify the calculation we can write : R = a + ib (6.58)

65)

where :

1 Z 2 0 + Z 4 1 tan(φ 1 ), and A 1 = 3(R 2 10 cos(2αd) + cos(2φ 1 )) + R 10 , B 2 = 2Z 1 Z 2 (R 2 10 cos(2φ 1 ) + cos(2φ 1 ) + R 10 ) + Z 2 2 , C 1 = 1 + R 2 10 + 2R 10 cos(2φ 1 ),

where K1 = 2Z 0 Z 1 tan(φ 1 ), K2 = 2Z 0 Z 3 1 tan(φ 1 ) , K3 = (Z 2 1 + Z 2 0 tan 2 (φ 1 )), K4 = (2Z 0 Z 2 1 + 2Z 0 Z 2 1 tan 2 (φ 1 )) , K5 = (Z 2 1 Z 2 0 + Z 4 1 tan 2 (φ 1 )) ; and K 1 = (sin(2αd) -R 2 10 sin(2φ 1 )), K 2 = Z 2 1 K 1 , K 3 = (R 2 10 + 2R 10 cos(2φ 1 ) + 1), K 4 = (2Z 1 -2Z 1 R 2 12 ), K 5 = Z 2 1 (R 2 10 -2R 10 cos(2φ 1 ) + 1). If we take this system : a -a = 0 b -b = 0 (6.68)

and by using this equality we find the following result :

The optimization method is the simplest to apply for this system to estimate the value of Z 2 , hence the expression of density :