Acoustic and geoacoustic inverse problems in randomly perturbed shallow-water environments

TL;DR

This paper develops a method to estimate the statistical properties of the water column and sea bottom in shallow-water environments using acoustic data, by modeling wave propagation as Markovian dynamics and solving a nonlinear inverse problem.

Contribution

It introduces a novel approach to estimate environmental parameters from acoustic data by formulating the problem as a Markovian and nonlinear inverse problem, validated with real data.

Findings

Successful estimation of water column fluctuations

Effective characterization of sea bottom properties

Validated method with experimental data

Abstract

The main goal of this paper is to estimate the regional acoustic and geoacoustic shallow-water environment from data collected by a vertical hydrophone array and transmitted by distant time-harmonic point sources. We aim at estimating the statistical properties of the random fluctuations of the index of refraction in the water column and the characteristics of the sea bottom. We first explain from first principles how acoustic wave propagation can be expressed as Markovian dynamics for the complex mode amplitudes of the sound pressure, which makes it possible to express the cross moments of the sound pressure in terms of the parameters to be estimated. We then show how the estimation problem can be formulated as a nonlinear inverse problem using this formulation, that can be solved by minimization of a misfit function. We apply this method to experimental data collected by the ALMA system (Acoustic Laboratory for Marine Applications).

Tables

Table 1. Table 1: Experimental and theoretical scintillation indices at different frequencies.

Equations

$$\displaystyle\Big{[}\rho(x,z)\nabla\cdot\Big{(}\frac{1}{\rho(x,z)}\nabla\Big{)}+\frac{\omega^{2}}{c(x,z)^{2}}\Big{]}\hat{p}(x,z)=\delta(x)\delta(z-z_{0}),$$

$$\displaystyle\Big{[}(\partial_{x}^{2}+\partial_{z}^{2})+\frac{\omega^{2}}{c(x,z)^{2}}\Big{]}\hat{p}(x,z)=\delta(x)\delta(z-z_{0}),$$

$$\Big{[}\partial_{r}^{2}+\frac{1}{r}\partial_{r}+\frac{1}{r^{2}}\partial_{\theta}^{2}+\partial_{z}^{2}+\frac{\omega^{2}}{c(r,z)^{2}}\Big{]}\hat{P}=\frac{1}{2\pi r}\delta(r)\delta(z-z_{0}).$$

$$\begin{array}[]{l}\bullet\mbox{ Dirichlet boundary condition on the top of the water column: }\mbox{$\quad\hat{p}(x,0)=0$ for all $x\in{\mathbb{R}}$,}\\ \bullet\mbox{ Continuity at depth $z_{\rm b}$:}\mbox{$\quad\hat{p}(x,z_{\rm b}^{-})=\hat{p}(x,z_{\rm b}^{+})$ and $\partial_{z}\hat{p}(x,z_{\rm b}^{-})/\rho_{\rm w}=\partial_{z}\hat{p}(x,z_{\rm b}^{+})/\rho_{\rm s}$.}\end{array}$$

$$\hat{p}{\bf 1}_{(0,+\infty)}(x)\in{\cal C}^{0}\big{(}(0,+\infty),H_{0}^{1}(0,+\infty)\cap H^{2}(0,+\infty)\big{)}\cap{\cal C}^{2}\big{(}(0,+\infty),L^{2}(0,+\infty)\big{)}.$$

$$c_0(z)=c_{\mathrm{w}}{1}_{[0,z_{\mathrm{b}}]}(z)+c_{\mathrm{s}}{1}_{(z_{\mathrm{b}},+\infty )}(z),$$

$$\mathcal{H}={\rho }_0(z){\partial }_z{\rho }_0(z{)}^{-1}{\partial }_z+\frac{{\omega }^2}{c_0(z{)}^2},$$

$$\begin{array}[]{l}\bullet\mbox{ Dirichlet at the top: $\phi(0)=0$,}\\ \bullet\mbox{ Continuity at $z=z_{\rm b}$: $\phi(z_{\rm b}^{+})=\phi(z_{\rm b}^{-})$ and $\phi^{\prime}(z_{\rm b}^{-})/\rho_{\rm w}=\phi^{\prime}(z_{\rm b}^{+})/\rho_{\rm s}$.}\end{array}$$

$$({\varphi }_1,{\varphi }_2)$$

$$={\rho }_{\mathrm{w}}^{-1}{\int }_0^{z_{\mathrm{b}}}\overline{{\varphi }_1(z)}{\varphi }_2(z)dz+{\rho }_{\mathrm{s}}^{-1}{\int }_{z_{\mathrm{b}}}^{+\infty }\overline{{\varphi }_1(z)}{\varphi }_2(z)dz.$$

$$(-\infty ,k_{\mathrm{s}}^2)\cup \{{\beta }_N^2,\dots ,{\beta }_1^2\},$$

$$\displaystyle\phi_{j}(z)=\left\{\begin{array}[]{l}A_{j}\sin(\sigma_{j}z/z_{\rm b})\mbox{ if }0\leq z\leq z_{\rm b}\\ A_{j}\sin(\sigma_{j})\exp(-\zeta_{j}(z-z_{\rm b})/z_{\rm b})\mbox{ if }z\geq z_{\rm b}\end{array}\right.$$

$${\sigma }_j=z_{\mathrm{b}}\sqrt{k_{\mathrm{w}}^2-{\beta }_j^2},{\zeta }_j=z_{\mathrm{b}}\sqrt{{\beta }_j^2-k_{\mathrm{s}}^2},$$

$$A_j^2=\frac{2/z_{\mathrm{b}}}{\frac{1}{{\rho }_{\mathrm{w}}}(1-\frac{\sin (2{\sigma }_j)}{2{\sigma }_j})+\frac{1}{{\rho }_{\mathrm{s}}}\frac{{\sin }^2({\sigma }_j)}{{\zeta }_j}}.$$

$$\tan (\sigma )=-\frac{\sigma }{\sqrt{z_{\mathrm{b}}^2(k_{\mathrm{w}}^2-k_{\mathrm{s}}{)}^2-{\sigma }^2}}\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{w}}},$$

$$\displaystyle\phi_{\gamma}(z)=\left\{\begin{array}[]{l}A_{\gamma}\sin(\eta_{\gamma}z/z_{\rm b})\mbox{ if }0\leq z\leq z_{\rm b}\\ A_{\gamma}\big{[}\sin(\eta_{\gamma})\cos(\xi_{\gamma}(z-z_{\rm b})/z_{\rm b})+\frac{\rho_{\rm s}}{\rho_{\rm w}}\frac{\eta_{\gamma}}{\xi_{\gamma}}\cos(\eta_{\gamma})\sin(\xi_{\gamma}(z-z_{\rm b})/z_{\rm b})\big{]}\mbox{ if }z\geq z_{\rm b}\end{array}\right.$$

$${\eta }_{\gamma }=z_{\mathrm{b}}\sqrt{k_{\mathrm{w}}^2-\gamma },{\xi }_{\gamma }=z_{\mathrm{b}}\sqrt{k_{\mathrm{s}}^2-\gamma }$$

$$A_{\gamma }^2=\frac{{\xi }_{\gamma }{\rho }_{\mathrm{s}}{z}_{\mathrm{b}}}{\pi ({\xi }_{\gamma }^2{\sin }^2({\eta }_{\gamma })+\frac{{\rho }_{\mathrm{s}}^2}{{\rho }_{\mathrm{w}}^2}{\eta }_{\gamma }^2{\cos }^2({\eta }_{\gamma }))}.$$

$$\left({\varphi }_{\gamma },\varphi \right)=\underset{M\to +\infty }{\lim }{\int }_0^M{\varphi }_{\gamma }(z)\varphi (z){\rho }_0(z{)}^{-1}dz,$$

$$\displaystyle\left(\phi,\phi\right)=\sum_{j=1}^{N}\big{|}\left(\phi_{j},\phi\right)\big{|}^{2}+\int_{-\infty}^{k_{\rm s}^{2}}\big{|}\left(\phi_{\gamma},\phi\right)\big{|}^{2}d\gamma.$$

$$\phi\mapsto\Big{(}\left(\phi_{j},\phi\right),j=1,\ldots,N,\left(\phi_{\gamma},\phi\right),\gamma\in({-\infty},{k_{\rm s}^{2}})\Big{)}$$

$$\Pi (r,r^{\prime })(\varphi )=\sum _{j=1}^N\left({\varphi }_j,\varphi \right){\varphi }_j{1}_{(r,r^{\prime })}({\beta }_j^2)+{\int }_r^{\min (k_{\mathrm{s}}^2,r^{\prime })}\left({\varphi }_{\gamma },\varphi \right){\varphi }_{\gamma }d\gamma 1_{(r,+\infty )}(k_{\mathrm{s}}^2),$$

$$\Pi (r,r^{\prime })(\mathcal{H}\varphi )=\sum _{j=1}^N{\beta }_j^2\left({\varphi }_j,\varphi \right){\varphi }_j{1}_{(r,r^{\prime })}({\beta }_j^2)+{\int }_r^{\min (k_{\mathrm{s}}^2,r^{\prime })}\gamma \left({\varphi }_{\gamma },\varphi \right){\varphi }_{\gamma }d\gamma 1_{(r,+\infty )}(k_{\mathrm{s}}^2).$$

$$\hat{p}(x,z)=\sum _{j=1}^N{\hat{p}}_j(x){\varphi }_j(z)+{\int }_0^{k_{\mathrm{s}}^2}{\hat{p}}_{\gamma }(x){\varphi }_{\gamma }(z)d\gamma +{\int }_{-\infty }^0{\hat{p}}_{\gamma }(x){\varphi }_{\gamma }(z)d\gamma .$$

$${\partial }_x^2{\hat{p}}_j+{\beta }_j^2{\hat{p}}_j$$

$${\partial }_x^2{\hat{p}}_{\gamma }+\gamma {\hat{p}}_{\gamma }$$

$$\hat{p}(x,z)=$$

$$+{\int }_{-\infty }^0\frac{{\hat{a}}_{\gamma ,0}}{|\gamma {|}^{1/4}}{e}^{-\sqrt{|\gamma |}x}{\varphi }_{\gamma }(z)d\gamma ,$$

$${\hat{a}}_{j,0}=$$

$${\hat{a}}_{\gamma ,0}=$$

$$c^2(x,z)=\frac{c_0^2(z)}{1+V(x,z)}.$$

$$V(x,z)=\nu (x,z)1_{(0,z_{\mathrm{b}})}(z)+i{\nu }_{\mathrm{s}}{1}_{(z_{\mathrm{b}},+\infty )}(z),$$

Full text

1 Introduction

In this paper we consider acoustic wave propagation in a randomly perturbed shallow-water waveguide with an absorbing sea bottom. The random perturbations of the index of refraction are due to internal waves, which induce temperature and salinity fluctuations, and the sea bottom is made of sediments, which induce dissipation. In the regime where random perturbations and dissipation are small and propagation distance is large, it is possible to get an effective description of the acoustic wave propagation in terms of a Markovian dynamics for the complex mode amplitudes of the expansion of the pressure field onto the guided modes of the unperturbed and non-dissipative waveguide. This Markovian dynamics involves coupling terms between guided modes and mode-dependent dispersion and loss terms. The coupling terms between guided modes come from the random perturbations in the water column. The effective dispersion and loss terms come from two effects: the random perturbations in the water column induce coupling between the guided and radiative modes and the deterministic dissipation in the sediment layer induce an exponential decay of the guided mode amplitudes. Both effects generate an irreversible, mode-dependent loss of energy carried by the guided modes.

The mathematical literature contains a lot of results on wave propagation in randomly perturbed waveguides motivated by underwater acoustics [18, 14, 16]. Those results derive from first principles and make it possible to relate the coefficients of the effective Markovian model to the physical parameters of the waveguide (in particular, the statistics of the fluctuations of the index of refraction and the complex acoustic impedance of the sea bottom). The mathematical statements also make it clear in which sense the Markovian model approximates the random dynamics of the mode amplitudes in the waveguide, because there are subtle effects that follow from the fact that the results are first established in a weak topology, as explained in Appendix A. Comparisons between theory and numerical simulations show good agreement, both for direct and synthetic inverse problems [2], but there is not so far much comparison between the detailed mathematical predictions and real experiments. On the other hand, the physical literature contains many theoretical results essentially based on coupled mode equations [1, 6, 7, 5, 3, 4]. However, the equations are derived ad hoc and it is not always straightforward to relate the coefficients of the effective equations to the physical and statistical parameters of the medium. This is problematic as we have in mind to use such results to solve an inverse problem.

As we will see in this paper, the effective attenuation of the mode amplitudes plays a key role. In the physical literature, modal attenuation coefficients are introduced directly in the coupled mode equations without derivation from first principles [8, 5, 17]. On the other hand most of mathematical studies do not take into account the attenuation of the bottom layer. However, even though dissipation is weak, it plays an important role in long-range propagation. For the direct problem point of view, the equipartition regime that is well characterized by Gaussian statistics [11] and a scintillation index (relative intensity variance) close to one in the absence of dissipation looks very different in the presence of weak attenuation and it may give rise to high fluctuations of the intensity. This was first pointed out by Creamer [5]. For the inverse problem point of view, signals measured on an hydrophone array can be processed to extract these attenuation coefficients, which in turn allow to get information on the medium, in particular, the bottom properties, as we will see in Section 5. We exploit data from at-sea experiments based on DGA’s ALMA (Acoustic Laboratory for Marine Applications) system [9]. We use data collected by a system consisting of four 32-element, vertical-line hydrophone array and a moored pinger located far away (at 9 km) and transmitting time-harmonic signals between $2$ and $13~{}kHz$ with a three-minute repetition rate. The measurements were carried out in November (7th-17th) 2016 near the shores of North East Corsica [10]. The signals recorded by the hydrophones can be processed by cross correlation calculations to estimate the properties of the medium. We explain this procedure and report its results in this paper.

The paper is organized as follows. The direct problem is analyzed in Sections 2-4. Section 3 is a review of the modal decomposition of the sound pressure in a homogeneous, non-dissipative waveguide. Section 4 describes the Markovian dynamics of the mode amplitudes in a random, dissipative waveguide. The inverse problem is formulated and solved using experimental data in Section 5.

2 Wave propagation in waveguides

Our model consists of a two-dimensional waveguide with range axis denoted by $x\in{\mathbb{R}}$ and transverse coordinate denoted by $z\in[0,+\infty)$. We suppose that the depth of the sea is constant and equals $z_{\rm b}$. When $z\in[0,z_{\rm b}]$ the medium is water, when $z>z_{\rm b}$ it becomes sediments. A point-like source at a fixed position $(x,z)=(0,z_{0})$ transmits a time-harmonic signal at frequency $\omega$ which is collected by a vertical array of receivers (hydrophones) at $x=x_{\rm a}$.

The Helmholtz equation for the acoustic pressure $\hat{p}(x,z)$ writes:

[TABLE]

for $x\in{\mathbb{R}},\;z\geq 0,$ where $\nabla=(\partial_{x},\partial_{z})$, $c(x,z)$, resp. $\rho(x,z)$, is the speed, resp. the density, at position $(x,z)$. We consider that $\rho$ is stepwise constant and equal to $\rho_{\rm w}$ in the water and $\rho_{\rm s}$ in the sediments. Therefore we have

[TABLE]

for $x\in{\mathbb{R}},\;z\in(0,z_{\rm b})\cup(z_{\rm b},+\infty)$.

Remark. The model (2) can be derived from a more realistic three-dimensional situation, in which the sound pressure field $\hat{P}$ satisfies in cylindrical coordinates

[TABLE]

The solution is radially symmetric and, if we neglect a near field factor of the form $\hat{p}/r^{5/2}$ then the scaled pressure field $\hat{p}(r,z)=\sqrt{r}\hat{P}(r,z)$ satisfies (2).

The differential equation (2) is endowed with the following boundary conditions:

[TABLE]

We are interested in solutions of (2) such that

[TABLE]

Let us now describe the results we obtain in this model. We compute in Section 3 the pressure for an ideal (homogeneous) waveguide in an exact form. We analyze how it behaves when the sound speed in the water is no longer constant but perturbed by some small noise and the sediments are weakly dissipative. The wave modes interact with each other and an asymptotic regime in the small noise and large distance limit is studied in Section 4. In particular we compute correlations of the recorded pressure signals for hydrophones located in a vertical segment for exponentially decaying correlations of the medium in Section 4.3.

3 Homogeneous waveguide

In this section, we consider a wave speed $c(x,z)$ which is constant both in the water and in the sediments:

[TABLE]

with $c_{\rm s}>c_{\rm w}$ (Pekeris waveguide). There is no dissipation, no fluctuation along the $x$-axis. The analysis of the perfect waveguide is classical [20] but we include it here for the sake of completeness. Let us introduce the Helmholtz operator

[TABLE]

where $\rho_{0}(z)=\rho_{\rm w}{\bf 1}_{[0,z_{\rm b}]}(z)+\rho_{\rm s}{\bf 1}_{(z_{\rm b},+\infty)}(z)$, with boundary conditions:

[TABLE]

We also denote $k_{0}(z)=k_{\rm w}{\bf 1}_{[0,z_{\rm b}]}(z)+k_{\rm s}{\bf 1}_{(z_{\rm b},+\infty)}(z)$, with $k_{\rm w}=\omega/c_{\rm w}$ and $k_{\rm s}=\omega/c_{\rm s}$. The Helmholtz operator is self-adjoint with respect to the scalar product defined on $L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$ by:

[TABLE]

The Helmholtz operator has a spectrum of the form

[TABLE]

where the $N$ modal wavenumbers $\beta_{j}$ are positive and $k_{\rm s}^{2}<\beta^{2}_{N}<\cdots<\beta_{1}^{2}<k_{\rm w}^{2}$.

Discrete spectrum. The $j$th eigenvector associated to the eigenvalue $\beta_{j}^{2}$ is

[TABLE]

where

[TABLE]

and

[TABLE]

The $\sigma_{j}$’s are the solutions in $(0,\sqrt{k_{\rm w}^{2}-k_{\rm s}^{2}}z_{\rm b})$ of

[TABLE]

and we denote by $N$ the number of solutions.

Continuous spectrum. For $\gamma\in(-\infty,k_{\rm s}^{2})$, the improper eigenvector has the form:

[TABLE]

where

[TABLE]

and

[TABLE]

We remark that $\phi_{\gamma}$ does not belong to $L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$, but $\left(\phi_{\gamma},\phi\right)$ can be defined for any test function $\phi\in L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$ as

[TABLE]

where the limit holds on $L^{2}(-\infty,k_{\rm s}^{2})$.

Completeness. We have for any $\phi\in L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$:

[TABLE]

The map which assigns to every element of $L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$ the coefficients of its spectral decomposition

[TABLE]

is an isometry from $L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$ onto ${\mathbb{C}}^{N}\times L^{2}(-\infty,k_{\rm s}^{2})$. There exists a resolution of the identity $\Pi$ such that for any $\phi\in L^{2}({\mathbb{R}}^{+},\rho_{0}^{-1})$ and $-\infty\leq r\leq r^{\prime}\leq+\infty$:

[TABLE]

and for any $\phi$ in the domain of ${\cal H}$,

[TABLE]

Modal decomposition. Any solution of the Helmholtz equation in homogeneous medium can be expanded as

[TABLE]

The first modes are guided, the second ones are radiating, the third ones are evanescent. Indeed, the complex mode amplitudes satisfy

[TABLE]

Therefore, if the source is in the plane $x=0$, at $(0,z_{0})$, as in (2), we have for $x>0$:

[TABLE]

with the modal amplitudes determined by the source:

[TABLE]

4 Random and dissipative waveguide

Here we assume that the waveguide is weakly randomly perturbed and weakly dissipative:

[TABLE]

The perturbation $V(x,z)$ has two components: random velocity perturbation in water and deterministic constant dissipation in sediments:

[TABLE]

where $\nu(x,z)$ is a zero-mean random process that describes the relative fluctuations of the propagation speed in water and $\nu_{\rm s}>0$ models the damping in the sediments.

4.1 The coupled mode equations

The solution of the perturbed Helmholtz equation

[TABLE]

for $x\in{\mathbb{R}},z\in(0,z_{\rm b})\cup(z_{\rm b},+\infty)$, with the boundary conditions (6), can be expanded as (22) and the complex mode amplitudes satisfy the coupled equations:

[TABLE]

with

[TABLE]

The coupling term $C_{jl}$ has two components:

[TABLE]

where

[TABLE]

and similarly for the other coupling terms.

We introduce the amplitudes of the generalized right and left-going mode amplitudes $\hat{a}_{j}$ and $\hat{b}_{j}$ for the guided and radiating modes:

[TABLE]

They satisfy the following coupled differential equations:

[TABLE]

[TABLE]

The radiating mode amplitudes $\hat{a}_{\gamma}(x)$ and $\hat{b}_{\gamma}(x)$, $\gamma\in(0,k_{\rm s}^{2})$, satisfy the same equations upon replacing $j$ by $\gamma$ in the formulas above. The evanescent mode amplitudes $\hat{p}_{\gamma}(x)$, $\gamma\in(-\infty,0)$, satisfy the differential equations

[TABLE]

where $\hat{g}_{\gamma}(x)$ and $\hat{g}_{\gamma}^{\rm ev}(x)$ are given by

[TABLE]

We wish to understand the evolution of the complex mode amplitudes in the regime where the perturbations are weak and the propagation distance is large. We assume

[TABLE]

where $\varepsilon$ is a small dimensionless parameter that expresses the fact that the random perturbations are weak and dissipation is even weaker. We suppose that we have a perfect (homogeneous and non-dissipative) waveguide outside the region $x\in(0,x_{\infty}/\varepsilon^{2})$ for some $x_{\infty}>0$. We carry out an asymptotic study in the small $\varepsilon$ limit. In this framework the evanescent mode amplitudes can be expressed in terms of the propagating mode amplitudes [13, 16] and the differential equation for $\hat{a}_{j}^{\varepsilon}(\omega,x):=\hat{a}_{j}(\omega,x/\varepsilon^{2})$ can be written as:

[TABLE]

In the asymptotic framework “weak fluctuations, weak dissipation, large propagation distance", we carry out a multiscale analysis similar as the one reported in Ref.  \onlinecitegomez. We find that the random fluctuations of the propagation speed in water give the standard coupling terms between the right-going mode amplitudes of the guided modes. They also induce loss towards the radiating modes. Moreover, the dissipation in the sediments gives rise to additional damping terms in the evolution equations of the complex mode amplitudes. The coupled system of equations for the complex mode amplitudes is of the form $d{\itbf a}^{\varepsilon}/dx=\varepsilon^{-1}{\itbf F}(x/\varepsilon^{2},{\itbf a}^{\varepsilon})+{\itbf G}(x/\varepsilon^{2},{\itbf a}^{\varepsilon})$. The small $\varepsilon$ limit of such differential equations are studied in the book  \onlineciteFGPSbook (Chapters 6 and 20). The full result for the complex mode amplitudes is expressed in Appendix A. From this result we get that the process $(|\hat{a}_{j}^{\varepsilon}(x)|^{2})_{j=1}^{N}$ converges towards a Markov process ${\itbf P}(x)=(P_{j}(x))_{j=1}^{N}$ whose infinitesimal generator $\mathcal{L}_{\itbf P}$ writes:

[TABLE]

with

[TABLE]

where

[TABLE]

and we have $\int_{z_{\rm b}}^{\infty}\phi_{j}(z)^{2}dz=A_{j}^{2}\sin^{2}(\sigma_{j})\,z_{\rm b}/(2\zeta_{j})$. This means that, for any test function $F:\mathbb{R}^{N}\to\mathbb{R}$, the expectation ${\mathbb{E}}[F({\itbf P}(x))]$ satisfies the Kolmogorov forward equation

[TABLE]

Equivalently, the probability density function $p(x,{\itbf P})$ of the random vector ${\itbf P}(x)$ satisfies the Kolmogorov backward (or Fokker-Planck) equation $\partial_{x}p={\cal L}_{\itbf P}^{*}p,$ where ${\cal L}_{{\itbf P}}^{*}$ is the adjoint of ${\cal L}_{\itbf P}$.

From the form of the generator ${\cal L}_{\itbf P}$, one can establish that the $n$th-order moments of the mode powers satisfy closed equations. In particular, using (45) and (50) we find that the mean mode powers

[TABLE]

satisfy the closed system of equations

[TABLE]

The form of these coupled-mode equations is well-known [6] although the mode-dependent attenuation term $\Lambda_{j}$ was usually introduced heuristically so far. The solution explicitly writes:

[TABLE]

with the matrix ${\bf A}$ defined by ($\delta_{jl}$ is the Kronecker symbol and $\Gamma_{jj}=-\sum_{l^{\prime}\neq j}\Gamma_{jl^{\prime}}$):

[TABLE]

4.2 Computation of the coefficients of the matrix ${\bf A}$

The goal of this paragraph is to show how to get closed-form expressions of the coefficients of the matrix ${\bf A}$ when the correlation of the perturbation decreases exponentially as a function of the horizontal distance between the two points, i.e.

[TABLE]

where $\ell_{\rm h}$ is the horizontal correlation radius of the random fluctuations of the index of refraction. We can compute $\Gamma_{jl}$ for $j\neq l$ from:

[TABLE]

and similarly we can compute $\Lambda_{j}$ from the same quantities upon substitution $\gamma$ for $l$. By (47) we have

[TABLE]

Replacing $\phi_{j}$ for all $j$ by their expression, and using the notation $k_{{\rm w}j}:=\sqrt{k_{\rm w}^{2}-\beta_{j}^{2}}$, we obtain:

[TABLE]

where

[TABLE]

When $R$ has the exponential form $R(z,z^{\prime})=\exp(-|z-z^{\prime}|/\ell_{\rm v})/2$, then the variance of the fluctuations is $\sigma^{2}/2$ and $S(k,k^{\prime})$ has a closed-form expression. We will use this particular form of the correlation function of the medium in Section 5 because it gives simple expressions for the coefficients $\Gamma_{jl}$ and $\Lambda_{j}$, which is convenient for the resolution of the inverse problem which requires many evaluations of such coefficients for different values $\sigma,\ell_{\rm v},\ell_{\rm h}$ of the statistics of the random medium and $\rho_{\rm s},\nu_{\rm s}$ of the sea bottom. We could as well use the Garrett-Munk correlation function at the expense of some computational overburden [15], but the impact of the exact form of the correlation function turns out to be weak.

4.3 Pressure field correlations

We now compute the correlation of the received signal on a vertical segment at a fixed horizontal distance $x_{\rm a}$. When $x_{\rm a}\gg 1$, we can use the asymptotic study from above and the decomposition of $\hat{p}(x_{\rm a},z)$ in the basis of eigenvectors $\phi_{j}$ is:

[TABLE]

Therefore

[TABLE]

The cross second moments ${\mathbb{E}}[\overline{\hat{a}_{j}}\hat{a}_{l}]$ for $j\neq l$ decay exponentially with the propagation distance $x$ as shown in Ref.  \onlineciteGP07 (see also Ref.  \onlinecitecolosi12) and they can be neglected as soon as the propagation distance $x$ is larger than the scattering mean free path, that depends on the coefficients $\Gamma_{jl}$. This is roughly speaking the distance beyond which the phase of the field has a random part which variance of the order of or larger than one, so that the coherent field is vanishing [14, 2]. The second moments ${\mathbb{E}}[|\hat{a}_{j}|^{2}]$ were computed in the previous paragraph using the exponential of the matrix ${\bf A}$:

[TABLE]

where $\hat{a}_{l,0}$ is given by $\eqref{a_j0}$. We can now write the correlation between the sound pressure signals recorded by two receivers which are separated by $y$ along the vertical segment at distance $x_{\rm a}$ and of depth $[z_{m},z_{M}]$ where $0<z_{m}<z_{M}<z_{\rm b}$. For all $y\in[0,z_{M}-z_{m}]$, the spatial correlation $C_{x_{\rm a}}(y)$ writes:

[TABLE]

4.4 Equipartition regime

By (53) the mean mode powers satisfy

[TABLE]

where $({\itbf V},-\lambda)$ is the first eigenvector/eigenvalue of the matrix ${\bf A}=\boldsymbol{\Gamma}-\boldsymbol{\Phi}$, with $\Gamma_{jl}$ given by (47) for $j\neq l$, $\Gamma_{jj}=-\sum_{l^{\prime}\neq j}\Gamma_{jl^{\prime}}$, $\Phi_{jl}=\Lambda_{j}\delta_{jl}$, and

[TABLE]

Note that:

• •

The coefficients of ${\itbf V}$ have all the same sign (so we can assume that they are nonnegative).

Indeed, the solution of the equation $\partial_{x}{\itbf Q}=(\boldsymbol{\Gamma}-\boldsymbol{\Phi}){\itbf Q}$ starting from any ${\itbf Q}(x=0)$ with nonnegative coefficients must have nonnegative coefficients since the exponential of a matrix whose off diagonal entries are nonnegative has nonnegative coefficients. As ${\itbf Q}(x)$ is equivalent to $({\itbf V}^{T}{\itbf Q}(x=0))\exp(-\lambda x){\itbf V}$, the coefficients of ${\itbf V}$ must have the same sign. We deduce as well that $c_{V}$ has the same sign as ${\itbf V}$.

• •

$\lambda\geq 0$.

Indeed, we have $(\boldsymbol{\Gamma}-\boldsymbol{\Phi}){\itbf V}=-\lambda{\itbf V}$. By projecting onto {\itbf V}^{(0)}=\big{(}1/\sqrt{N}\big{)}_{j=1}^{N}, which is in the kernel of $\boldsymbol{\Gamma}$, we get ${{\itbf V}^{(0)}}^{T}\boldsymbol{\Phi}{\itbf V}=\lambda{{\itbf V}^{(0)}}^{T}{\itbf V}$. Since $\boldsymbol{\Phi}$ is diagonal with nonnegative coefficients and the coefficients of ${\itbf V}$ have the same sign and cannot be all zero, we find $\lambda=\sum_{j}\Lambda_{j}{V}_{j}/\sum_{j}{V}_{j}\geq 0$.

In the following we discuss cases with zero or weak dissipation, where explicit expressions can be obtained. Remember that the effective dissipation is the sum of the two effects: power leakage from the guided modes to the radiating modes and dissipation in the sediments.

No effective dissipation. If there is no effective dissipation $\boldsymbol{\Phi}={\bf 0}$, then the first eigenvector/eigenvalue $({\itbf V}^{(0)},-\lambda^{(0)})$ of the matrix $\boldsymbol{\Gamma}$ is

[TABLE]

which gives the standard equipartition result [5, 12, 14]:

[TABLE]

Weak effective dissipation. We next consider the case when the effective dissipation is weak, that is to say, the matrix $\boldsymbol{\Phi}$ is much smaller than the matrix $\boldsymbol{\Gamma}$, with a typical ratio of the order of $\delta$. We then assume that $\Lambda_{j}=\delta\Lambda_{j}^{(1)},$ with $\delta\ll 1$. Then we can write $\boldsymbol{\Phi}=\delta\boldsymbol{\Phi}^{(1)}$ with $\Phi_{jl}^{(1)}=\Lambda_{j}^{(1)}\delta_{jl}$ and $\boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{(0)}$ and the first eigenvector/eigenvalue $({\itbf V},-\lambda)$ of the matrix $\boldsymbol{\Gamma}-\boldsymbol{\Phi}$ can be expanded as

[TABLE]

with

[TABLE]

and ${\itbf V}^{(1)}$ is solution of $\boldsymbol{\Gamma}^{(0)}{\itbf V}^{(1)}=(\boldsymbol{\Phi}^{(1)}-\lambda^{(1)}){\itbf V}^{(0)}$ and is orthogonal to ${\itbf V}^{(0)}$. If, for instance, $\Gamma_{jl}\equiv\Gamma>0$ for all $j\neq l$, then

[TABLE]

and

[TABLE]

4.5 Fluctuation analysis

By (45) we find that the second-order moments of the mode powers

[TABLE]

satisfy the closed equations

[TABLE]

for $j\neq l$. This system has the same form as the one found in the literature dedicated to coupled mode theory [6, 5]. The initial conditions are $R_{jl}(0)=|\hat{a}_{j,0}|^{2}|\hat{a}_{l,0}|^{2}$. Let us introduce ${\itbf S}=(S_{jl})_{1\leq j\leq l\leq N}$ defined by

[TABLE]

The $S_{jl}$’s satisfy the system

[TABLE]

with the convention that whenever $S_{jl}$ occurs with $j>l$, it is replaced by $S_{lj}$. This can be written in the form $\partial_{x}{\itbf S}=(\boldsymbol{\Theta}-\boldsymbol{\Psi}){\itbf S}$. The linear operator $\boldsymbol{\Psi}$ (that depends on the $\Lambda_{j}$’s) is diagonal and the linear operator $\boldsymbol{\Theta}$ (that depends on the $\Gamma_{jl}$’s) is self-adjoint: for any ${\itbf T}$ and $\widetilde{\itbf T}$, we have

[TABLE]

because $2\sum_{j<l}=\sum_{j\neq l}$. As a consequence, $\boldsymbol{\Theta}-\boldsymbol{\Psi}$ can be diagonalized and we find that

[TABLE]

where $c_{W}$ is the projection on the first eigenvector ${\itbf W}$ in the basis of eigenvectors of $\boldsymbol{\Theta}-\boldsymbol{\Psi}$, and

[TABLE]

with the convention that whenever $W_{jl}$ occurs with $j>l$, it is replaced by $W_{lj}$.

No effective dissipation. If there is no effective dissipation, then the first eigenvector/eigenvalue $({\itbf W}^{(0)},-\mu^{(0)})$ of the matrix $\boldsymbol{\Theta}$ is

[TABLE]

with $c_{N}=\sqrt{2}/\sqrt{N(N+1)}$ and

[TABLE]

As $\sum_{j\leq l}S_{jl}(x)=\sum_{j,l}R_{jl}(x)=(\sum_{j=1}^{N}|\hat{a}_{j,0}|^{2})^{2}$, we deduce

[TABLE]

and

[TABLE]

Weak effective dissipation. We next consider the case when the effective dissipation is weak, say $\Lambda_{j}=\delta\Lambda_{j}^{(1)}$. Then we can write $\boldsymbol{\Psi}=\delta\boldsymbol{\Psi}^{(1)}$ and $\boldsymbol{\Theta}=\boldsymbol{\Theta}^{(0)}$ and the first eigenvector/eigenvalue $({\itbf W},-\mu)$ of the matrix $\boldsymbol{\Theta}-\boldsymbol{\Psi}$ can be expanded as

[TABLE]

with

[TABLE]

and ${\itbf W}^{(1)}$ is solution of $\boldsymbol{\Theta}^{(0)}{\itbf W}^{(1)}=(\boldsymbol{\Psi}^{(1)}-\mu^{(1)}){\itbf W}^{(0)}$ and is orthogonal to ${\itbf W}^{(0)}$. If, for instance, $\Gamma_{jl}\equiv\Gamma>0$ for all $j\neq l$, then

[TABLE]

and

[TABLE]

Note that

[TABLE]

is negative-valued.

Exponential growth of the intensity fluctuations. It is a general feature that, for any matrix $\boldsymbol{\Gamma}$ and effective dissipation coefficients $\Lambda_{j}$, we have $\mu-2\lambda\leq 0$ (we have equality when there is no effective dissipation and this is a consequence of the forthcoming result (67) and Cauchy-Schwarz inequality when there is dissipation). The first two moments of the pointwise intensity $|\hat{p}(x,z)|^{2}$ for large $x$ are

[TABLE]

Without dissipation we have the following result for the fluctuations of the pointwise intensity:

[TABLE]

which is equal to $2$ when $N\gg 1$, and with dissipation

[TABLE]

that grows exponentially with the propagation distance (for very long distances, however, as $\mu-2\lambda$ is very small as shown above). Eq. (64) gives the expression of the exponential growth rate when dissipation is weak. In this regime the growth rate increases when the effective modal dissipation coefficients become different from each other and decreases when the number of modes increases.

5 Inverse problem

The goal of this section is to show that it is possible to estimate the statistical properties of the index of refraction in the water column and the sea bottom properties from the sound pressure recorded by a vertical hydrophone array and transmitted by distant time-harmonic sources. This inverse problem can be formulated as a minimization problem that tries to match empirical quantities with theoretical ones that depend on the parameters to be estimated. The theoretical model of the previous section shows that the correlation function of the sound pressure has a non-trivial behavior that makes it possible to identify many relevant parameters as we explain below. We first present the ALMA experimental data, then formulate the inverse problem, and finally estimate the model parameters.

5.1 ALMA 2016 experiment

ALMA is a series of at-sea experiments carried out by French DGA Naval Systems [9]. In 2016, the experiment took place in november 7th-17th near the shores of North East Corsica [10] (cf. Figure 1). There was a pinger immerged at fifty meters, and nine kilometers away, a passive array with $128$ hydrophones ($4$ verticals arms of $32$ hydrophones, the arms are $0.5~{}m$ apart from each other, the hydrophones are $0.15~{}m$ apart from each other) immerged at sixty meters. The sea bottom is relatively flat on the line between the pinger and the array (between $100~{}m$ and $115~{}m$), and composed of sand and gravelly-sand. Sea was calm, with roughness height around $0.1~{}m$ (sea state 1). The pinger transmitted, for about 2 hours and a half, a sequence consisting of several time-harmonic waves, with a repetition rate of three minutes. The sequence was a train made of two-second time-harmonic waves at the following $K=6$ frequencies: $2kHz$, $5~{}kHz$, $7~{}kHz$, $9~{}kHz$, $11~{}kHz$, $13~{}kHz$. Acquisition was sampled at $48~{}kHz$. The Fourier coefficients for each frequency were extracted using a Hann window function with a duration of $1~{}s$.

Each arm has the form of a vertical array at a distance $x_{\rm a}$ of the source immerged at $z_{0}$. The hydrophones are regularly spaced between the depths $z_{m}$ and $z_{M}$, with $0<z_{m}<z_{M}<z_{\rm b}$. At frequency $\omega$, for $y\in\left[0,z_{M}-z_{m}\right]$, we denote the spatial correlation between hydrophones at a distance $y$ from each other by $C_{x_{\rm a}}(y)$ (cf. Eq. (55) for the theoretical expression predicted by our model). Then we define the correlation radius $r$ as the half-width at half-maximum (i.e. the distance between hydrophones for which the correlation is $1/2$):

[TABLE]

From the experimental data at frequencies $f_{1},\dots,f_{K}$ (with $\omega=2\pi f$) we can extract the experimental correlation radii $r_{\rm e}(f_{1}),\dots,r_{\rm e}(f_{K})$. For $\Phi$ a set of parameters of our theoretical model (assuming that we know the sound speed $c_{\rm w}$ and density $\rho_{\rm w}$ in water and the depth of the sea bottom),

[TABLE]

we can compute the theoretical correlation radii $r_{\rm t}(f_{1},\Phi),\dots,r_{\rm t}(f_{K},\Phi)$. Then we can define the misfit function

[TABLE]

and we can determine the parameters $\Phi$ by minimizing over $\Phi$ the misfit function $E(\Phi)$.

5.2 Estimation of the model parameters

We assume that sound speed is constant in water, taking value $c_{\rm w}=1523m/s$ from CTD (conductivity, temperature, and depth) measurements, and we set density of water $\rho_{\rm w}$ at $1000~{}kg/m^{3}$ (we could have taken a more precise value, but density has little influence on the model). Here, we have $x_{\rm a}=9000~{}m$. Both source and array are in the middle of the water column. Figure 2 presents the correlation radii obtained from experimental data for each frequency, and the theoretical correlation radii predicted by our model with the parameters determined by minimization of the misfit function (68): $c_{\rm s}=1630~{}m/s$; $\alpha=1.09~{}dB/\textrm{wavelength}$; $\sigma=0.002$ (relative fluctuation); $\ell_{\rm v}=30~{}m$; $\ell_{\rm h}=100~{}m$; $\rho_{\rm s}=1700~{}kg/m^{3}$. Parameter $\nu_{\rm s}$ is computed from the attenuation coefficient $\alpha$ expressed in decibel by wavelength. The values of the parameters seem compatible with experimental measurements carried out in similar environments [19]. Figure 3 shows theoretical and experimental correlation functions at the different frequencies.

We carry out a basic (one-at-a-time) sensitivity analysis: we modify one by one the optimal values we have found, letting the others unchanged. Results are in Figure 4. We see that $c_{\rm s}$ has a noticable influence, as a $20~{}m/s$ change in sediment sound speed significantly modifies the correlation radii. $\sigma$ is also important, as dividing or multiplying it by two has a strong influence. $\alpha$ and $\ell_{\rm v}$ have less influence, but are still significant. Finally, for reasonable physical changes of $\rho_{\rm s}$ and $\ell_{\rm h}$, we do not see noticeable effects in the results, whatever the frequency is.

The inspection of the numerical values of the two terms that determine $\Lambda_{j}$ (cf. Eq. (4.1)) give interesting information. We can see that for each frequency, the radiative terms are approximately ten times smaller than the attenuation ones. So we can claim that, in the experimental configuration addressed in this paper, radiation effects can be neglected.

Finally, if we denote by $(z_{n})_{1\leqslant n\leqslant N_{\rm h}}$ the depths of the $N_{\rm h}$ hydrophones of the array, we can compare the theoretical frequency-dependent scintillation indices defined by

[TABLE]

with the experimental values $S_{\rm e}$ determined by empirical averages instead of expectations. Table 1 shows experimental scintillation indices compared to theoretical ones, with the optimal parameters $\Phi$ determined from the correlation radii. In the experimental and theoretical configurations scintillation indices are approximately equal to one. One should have longer propagation distances to observe scintillation indices larger than one, as already suggested in Ref.  \onlinecitecreamer.

6 Conclusion

This paper proposes a complete description of the statistics of the mode amplitudes of the sound pressure in a shallow-water waveguide. The effective parameters (frequency- and mode-dependent attenuation, dispersion, and coupling) are identified from first principles and expressed in terms of the statistical properties of the index of refraction of the water column and the sea bottom properties. This theoretical analysis makes it possible to formulate an inverse problem for the estimation of the model parameters from the sound pressure recorded by a vertical hydrophone array and transmitted by distant time-harmonic sources. This inverse problem is solved using data collected during the ALMA 2016 experiment.

In the experimental configuration addressed in this paper, it turns out that radiation effects can be neglected (compared to dissipation in the sediments) and non-Gaussian scintillation effects can be neglected as well. The extraction of the frequency-dependent correlation radii of the recorded sound pressure signals makes it possible to estimate different acoustic and geoacoustic parameters: The sediment sound speed $c_{\rm s}$ and the standard deviation of the index of refraction in the water column $\sigma$ can be robustly estimated because they have strong effects on the correlation radii of the sound pressure. The attenuation parameter in the sediments $\alpha$ and the vertical correlation radius of the index of refraction in the water column $\ell_{\rm v}$ can be estimated, because they have less influence, but are still significant. The sediment density $\rho_{\rm s}$ and the horizontal correlation radius $\ell_{\rm h}$ are difficult to estimate because they have small effects of the correlation radii of the sound pressure.

Acknowledgements

We thank Dominique Fattaccioli and Gaultier Real from DGA Naval Systems for interesting and stimulating discussions and for making us available the ALMA 2016 data.

Appendix A The effective Markovian system for the complex mode amplitudes

We assume that the reduced wavenumbers $\beta_{j}$ are distinct, for $1\leq j\leq N$. Then, for any $x_{\infty}>0$, the process \big{(}(\hat{a}_{j}^{\varepsilon}(x))_{j=1}^{N},(\hat{a}_{\gamma}^{\varepsilon}(x))_{\gamma\in(0,k_{\rm s}^{2})}\big{)} converges in distribution in ${\cal C}^{0}([0,x_{\infty}],\mathbb{C}^{N}\times L^{2}(0,k_{\rm s}^{2}))$, where $\mathbb{C}^{N}\times L^{2}(0,k_{\rm s}^{2})$ is equipped with the weak topology, to the Markov process \big{(}({a}_{j}(x))_{j=1}^{N},({a}_{\gamma}(x))_{\gamma\in(0,k_{\rm s}^{2})}\big{)} with infinitesimal generator ${\cal L}={\cal L}^{1}+{\cal L}^{2}+{\cal L}^{3},$ where ${\cal L}^{j}$, $1\leq j\leq 3$, are the differential operators:

[TABLE]

In these definitions we use the classical complex derivative: if $\zeta=\zeta_{r}+i\zeta_{i}$, then $\partial_{\zeta}=(1/2)(\partial_{\zeta_{r}}-i\partial_{\zeta_{i}})$ and $\partial_{\overline{\zeta}}=(1/2)(\partial_{\zeta_{r}}+i\partial_{\zeta_{i}})$, and the coefficients of the operators (69-71) are defined for indices $j,l=1,\ldots,N$, as follows:

For all $j\neq l$, $\Gamma_{jl}$ is defined by (47) and

[TABLE]

For all $j,l$:

[TABLE]

For all $j$, $\Lambda_{j}$ is defined by (4.1) and

[TABLE]

Remarks:

1. The convergence result holds in the weak topology. This means that we can only compute quantities of the form ${\mathbb{E}}[F(a_{1},\ldots,a_{N},\int_{0}^{k_{\rm s}^{2}}\alpha_{\gamma}a_{\gamma}d\gamma)]$ for any test function $\alpha\in L^{2}(0,k_{\rm s}^{2})$ and $F:\mathbb{R}^{N+1}\to\mathbb{R}$, which are the limits of ${\mathbb{E}}[F(\hat{a}^{\varepsilon}_{1},\ldots,\hat{a}^{\varepsilon}_{N},\int_{0}^{k_{\rm s}^{2}}\alpha_{\gamma}\hat{a}^{\varepsilon}_{\gamma}d\gamma)]$ as $\varepsilon\to 0$.

2. The generator ${\cal L}$ does not involve $\partial_{a_{\gamma}}$. Therefore $(\hat{a}_{j}^{\varepsilon}(x))_{j=1}^{N}$ converges in distribution in ${\cal C}^{0}([0,x_{\infty}],\mathbb{C}^{N})$ to the Markov process $({a}_{j}(x))_{j=1}^{N}$ with generator ${\cal L}$. The weak and strong topologies are the same in $\mathbb{C}^{N}$, so we can compute any moment of the form ${\mathbb{E}}[F(a_{1},\ldots,a_{N})]$, which are the limits of ${\mathbb{E}}[F(\hat{a}^{\varepsilon}_{1},\ldots,\hat{a}^{\varepsilon}_{N})]$.

3. ${\cal L}_{1}$ is the contribution of the coupling between guided modes, which gives rise to power exchange between the guided modes; ${\cal L}_{2}$ is the contribution of the coupling between guided and radiating modes, which gives rise to power leakage from the guided modes to the radiating ones (effective diffusion) and addition of frequency-dependent phases on the guided modes (effective dispersion); ${\cal L}_{2}$ also contains the effective mode-dependent term due to dissipation in the sediments; ${\cal L}_{3}$ is the contribution of the coupling between guided and evanescent modes, which gives rise to additional phase terms on the guided modes (effective dispersion [13]).

4. If the generator ${\cal L}$ is applied to a test function that depends only on the mode powers $(P_{j}=|{a}_{j}|^{2})_{j=1}^{N}$, then the result is a function that depends only on $(P_{j})_{j=1}^{N}$. Thus, the mode powers $(P_{j}(x))_{j=1}^{N}$ define a Markov process, with infinitesimal generator defined by (45).

5. The radiation mode amplitudes remain constant on $L^{2}(0,k_{\rm s}^{2})$, equipped with the weak topology, as $\varepsilon\to 0$. However, this does not describe the power $\int_{0}^{k_{\rm s}^{2}}|\hat{a}^{\varepsilon}_{\gamma}|^{2}d\gamma$ transported by the radiation modes, because the convergence does not hold in the strong topology of $L^{2}(0,k_{\rm s}^{2})$ so we do not have $\int_{0}^{k_{\rm s}^{2}}|\hat{a}^{\varepsilon}_{\gamma}|^{2}d\gamma\to\int_{0}^{k_{\rm s}^{2}}|a_{\gamma}|^{2}d\gamma$ as $\varepsilon\to 0$.