<figure><img src="image_1.png" /></figure>

HAL
open science

# Sound-Source Localization in Range-Dependent Shallow-Water Environments Using a Four-Layer Model

Xun Wang, Shahram Khazaie, Dimitri Komatitsch, Pierre Sagaut

► To cite this version:

Xun Wang, Shahram Khazaie, Dimitri Komatitsch, Pierre Sagaut. Sound-Source Localization in Range-Dependent Shallow-Water Environments Using a Four-Layer Model. IEEE Journal of Oceanic Engineering, 2017, 44 (1), pp.220 - 228. 10.1109/JOE.2017.2775978. hal-01702364

HAL Id: hal-01702364

https://hal.science/hal-01702364v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sound-Source Localization in Range-Dependent
Shallow-Water Environments Using
a Four-Layer Model

Xun Wang, Shahram Khazaie, Dimitri Komatitsch, and Pierre Sagaut

*Abstract*—Sound-source localization in shallow water is a diffi-cult task due to the complicated environment, e.g., complex sound-speed profile and irregular water bottom reflections. Full-wave numerical techniques are currently able to accurately simulate the propagation of sound waves in such complex environments. How-ever, the source localization problem, which generally involves a large number of sound propagation calculations, still requires a fast computation of the wave equation, and thus a simplified model is well advised. In this paper, a four-layer model is considered, which is able to approximate a wide range of shallow-water envi-ronments, particularly those in summer conditions. More specif-ically, the medium is assumed to be horizontally stratified and vertically divided into four layers, and the sound speed in each layer is assumed to be constant or varying linearly. Under this assumption, the wave propagation can be rapidly computed via a classical wave number integration method. The main contribution of this paper is to show the suitability of the four-layer model in terms of source localization in a complex (range-dependent) envi-ronment. The sound-speed profile is assumed to be vertically irreg-ular and horizontally slowly varying and the bottom is nonflat. In the forward problem, sound propagation in complex underwater environments is simulated via a time-domain full-wave simulation approach called the spectral-element method. The source localiza-tion error due to model imprecision is analyzed.

*Index Terms*—Matched-field processing, range-dependent environment, shallow water, sound-source localization, spectral-element method.

## I. INTRODUCTION

SOURCE localization in shallow water is a chal-lenging problem due to the complexity of the underwater environment. The complex nature of the sound-speed distribu-tion, which may be heterogeneous and imprecisely known, is

especially important [1]–[8]. Particularly, source localization in range-dependent environments [9]–[14], which may involve ocean bottom slopes and horizontal variations in sound speed and density due to, for example, salinity, temperature, and internal waves, has become of significant interest in recent years.

To precisely localize a sound source, the first step consists in modeling the sound propagation in the underwater envi-ronment. For this purpose, various numerical methods can be employed; for example, finite-difference methods, boundary-element methods, classical low-order or spectral-element meth-ods (SEMs) have proven to be efficient tools to simulate sound propagation in a realistic and complicated medium [15], [16]. However, for a high-frequency source, to be accurate, numerical methods require a large number of mesh elements. Furthermore, in the underwater source localization problem, the sound prop-agation range is long and signal processing techniques require a large number of computations of the sound field. The above reasons imply that the numerical methods involve high com-putational costs. Other solutions for simulating wave propaga-tion in a range-dependent environment are based on the normal mode approach [17] or parabolic equations [9], [10], [18]. How- ever, these techniques still involve a large amount of numerical computations and are constrained to some specific cases; for example, the environment can be divided into a sequence of range-independent regions. For the above reasons, the model used for sound field computation is always simplified; for ex-ample, using a planar multilayered model in which the sound speed varies linearly in each layer (range independent) [19]–[23], such that the wave propagation can be rapidly computed based on simpler techniques.

In this paper, a four-layer model is considered for model ap-proximation in the inverse problem by its availability in terms of giving accurate source localization in complex range-dependent environments is studied. The sound speed is assumed to be hor-izontally stratified; the vertical profile in the water column con-sists of two separate constant layers, between which the sound speed linearly varies. Moreover, a semi-infinite fluid bottom is considered, in which the sound speed is constant. The densi-ties in the water column and ocean bottom are both constant but with different values. It has been mentioned in [1] that this sound-speed profile is able to approximate a kind of shallow water medium, particularly those during the summer months. The depth-separated wave equation under this assumption can be analytically solved [19]. A rapid computation of sound

X. Wang is with the Aix-Marseille University, CNRS, Centrale Marseille, M2P2, Marseille 13451, France and also with the Department of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Hong Kong (e-mail: xunwang00@gmail.com).

S. Khazaie is with the Aix-Marseille University, CNRS, Centrale Mar-seille, M2P2, Marseille 13451, France and also with the Institut de Recherche Technologique Jules Verne, Chemin du Chaffault, Bouguenais 44340, France (e-mail: shahramkhazaee@gmail.com).

D. Komatitsch is with the Aix-Marseille University, CNRS, Centrale Mar-seille, LMA, Marseille 13453, France (e-mail: komatitsch@lma.cnrs-mrs.fr).

P. Sagaut is with the Aix-Marseille University, CNRS, Centrale Marseille, M2P2, Marseille 13451, France (e-mail: pierre.sagaut@univ-amu.fr).
propagation between any two points in the medium can thus be obtained, which only needs to calculate a onefold numerical integration (inverse Hankel transform) based on the wave number integration method [15], [19], [20].

This paper studies the applicability of this simplified model in terms of source localization in a realistic environment. More specifically, the localization error in a range-dependent environment is quantified, e.g., nonflat ocean bottom or horizontally varying sound speed. In the forward problem, the sound propagation in a range-dependent environment is simulated via a time-domain full-wave simulation technique called the SEM [24], [25]. This approach is a formulation of the finite-element method (FEM) that uses high-degree polynomials as elemental basis functions. The accuracy of this numerical method for sound wave simulation in the ocean has been validated [16]. The underwater environment is considered in a 3-D space but the wave equation is assumed to be independent of the horizontal azimuth, thus the 2-D axisymmetric version of SEM can be used [26]. We thus use the open-source code SPECFEM2D [16], which is available at geodynamics.org, to simulate the sound wave propagation in the time domain.

In this paper, two matched-field processing (MFP) approaches [2]–[5], [15], [27]–[30] are used as the signal processing techniques for source localization. Conventional MFP [31], [32] maximizes the power output of a point source, which leads to a maximum likelihood estimate of the source location and strength [33]–[35]. Minimum variance distortionless response filter, or Capon's MFP [36], [37], minimizes the variance at the output of a linear weighting of the sensors subject to unit gain. Compared to the conventional MFP, this approach returns a superresolution estimate and compresses the level of side lobes, but it needs a relatively large number of samples to guarantee the full rank of sample covariance matrix and is more sensitive to environmental parameter uncertainties. More detailed introduction for the two MFP approaches and other versions of MFP, as well as their comparison can be found in review papers [27], [33].

The plan of this paper is as follows. Section II introduces the four-layer model and the solution of wave equation for this model. The sound-source localization methods are briefly presented in Section III. Then, numerical examples are introduced in Section IV. A summer Mediterranean shallow-water sound speed is approximated by the four-layer model, and the influence of nonflat bottom and horizontally varied sound speed on the source localization is investigated. Finally, conclusions are drawn in Section V.

## II. SOUND PROPAGATION IN A FOUR-LAYER SHALLOW WATER ENVIRONMENT

In this section, the sound propagation in shallow water with a four-layer model is considered, as illustrated in Fig. 1. A sound source $\mathbf{r}_0 = (x_0, y_0, z_0)$ is placed in the underwater region $\{\mathbf{r} = (x, y, z) : z \in [0, d)\}$, where $d > 0$ stands for the water depth (the water bottom is assumed to be flat for now in this simplified model). The sound speed in this region is horizontally stratified, continuous with respect to $z$, and in the vertical direction consists of two separate constant layers, with a layer in between in which

<figure><img src="image_3.png" /><figcaption>Fig. 1. Shallow-water environment with a four-layer sound-speed profile $c(z)$. The cross and circles represent the source and measurement locations, respectively.</figcaption></figure>

$k^2(z)$ varies linearly with respect to depth (the sound speed $c(z)$ is thus pseudolinear [15]). Here, $k(z) = 2\pi f/c(z)$ is the wave number and $f$ is the frequency. The bottom is represented by an infinite fluid half-space $\{\mathbf{r} = (x, y, z) : z \in [d, \infty)\}$ in which the sound speed is constant. The sound-speed profile is shown on the left of Fig. 1 and is represented as follows:

$$c(z) = \begin{cases} c_1, & z \in [0, d_1) \\ \sqrt{\frac{1}{az + b}}, & a > 0, z \in [d_1, d_2) \\ c_2, & z \in [d_2, d) \\ c_3, & z \in [d, \infty) \end{cases} \quad (1)$$

where the coefficients $a$ and $b$ are

$$a = \frac{c_1^2 - c_2^2}{c_1^2 c_2^2 (d_2 - d_1)}, \quad b = \frac{c_2^2 d_2 - c_1^2 d_1}{c_1^2 c_2^2 (d_2 - d_1)}. \quad (2)$$

In both the water column and the ocean bottom, the density is assumed to be constant

$$\rho(z) = \begin{cases} \rho_1, & z \in [0, d) \\ \rho_2, & z \in [d, \infty). \end{cases} \quad (3)$$

The wave equation for the displacement potential, denoted by $\tilde{\psi}$, as a function of the spatial coordinate $\mathbf{r} = (x, y, z)$ and time $t$, which exists because the density is constant in each layer [15], is governed by

$$\left(\nabla^2 - \frac{1}{c^2(z)} \frac{\partial^2}{\partial t^2}\right) \tilde{\psi}(\mathbf{r}, t) = \delta(\mathbf{r} - \mathbf{r}_0)S_t \quad (4)$$

where $S_t$ is the deterministic function of the source signal in the time domain, and $\delta$ is the Dirac delta distribution. Taking a Fourier transform of both sides of (4) with respect to $t$ results in the Helmholtz equation for the sound field $\psi(\mathbf{r}, f)$ in the frequency domain

$$(\nabla^2 + k^2(z)) \psi(\mathbf{r}, f) = \delta(\mathbf{r} - \mathbf{r}_0) S_f \quad (5)$$

where $S_f$ is the Fourier transform of $S_t$. The sound pressure is obtained from the displacement potential as
$$p(\mathbf{r}, f) = \rho\omega^2\psi(\mathbf{r}, f) \qquad (6)$$

where $\omega = 2\pi f$ is the angular frequency.

Note that by assuming that the sound source is omnidirectional, the sound field produced by the source only depends on the depth and on the horizontal range. For this reason, a cylindrical coordinate system is chosen, i.e., the spatial coordinate is denoted by $\mathbf{r} = (r, z, \phi)$, where $r$ and $\phi$ stand for the horizontal range and azimuth, respectively. By applying the Hankel transform

$$f(k_r, z) = \int_0^\infty f(r, z) J_0(k_r r) r dr \qquad (7)$$

to (5), the depth-separated wave equation is obtained

$$\left[ \frac{\partial^2}{\partial z^2} + (k^2(z) - k_r^2) \right] \psi(k_r, z) = S_f \frac{\delta(z - z_0)}{2\pi} \qquad (8)$$

where $k_r$ is the horizontal wave number. Note that (8) can be solved analytically, as recalled in Appendix A. Then, by applying an inverse Hankel transform on $\psi(k_r, z)$, i.e.,

$$\psi(r, z) = \int_{0}^{\infty} \psi(k_{r}, z) J_{0}(k_{r} r) k_{r} dk_{r} \qquad (9)$$

the sound field at any point can be computed. Here, the integral (9) is computed numerically based on the wave number integration method as recalled in Appendix B.

### III. SOUND-SOURCE LOCALIZATION

In this section, the MFP methods for solving the source localization problem are briefly recalled. As shown in Fig. 1, $M$ measurement stations are placed at $\mathbf{r}'_m = (x'_m, y'_m, z'_m)$, $m = 1, \dots, M$, and the sound pressure measurements are denoted by $\mathbf{p}_t = (p_{1t}, \dots, p_{Mt})^\top$, where $t = 1, \dots, T$ is the number of snapshots and $\top$ stands for the transpose operator. The Green's function for pressure (with unit amplitude), which corresponds to the sound propagation from the source $\mathbf{r}_0$ to $\mathbf{r}$, is obtained from (5) and (6). Let $G(\mathbf{r}, \mathbf{r}_0)$ denote the Green's function and $G(\mathbf{r}_0) = (G(\mathbf{r}_1, \mathbf{r}_0), \dots, G(\mathbf{r}_M, \mathbf{r}_0))^\top$ represent the vector of Green's functions at the measurement locations. Then, the acoustical measurement has the theoretical expression

$$\mathbf{p}_t = S\mathbf{G} + \mathbf{n}_t \qquad (10)$$

where $S$ is the amplitude of the source signal and $\mathbf{n}_t$ is an $M$-dimensional complex-valued random vector representing the measurement error. The sample covariance matrix can be estimated based on the $T$ snapshots: $\hat{\mathbf{K}} = \frac{1}{T} \sum_{t=1}^T \mathbf{p}_t \mathbf{p}_t^\mathrm{H}$ where $H$ is the Hermitian transpose. The sound source can then be localized using the conventional and Capon's MFP methods [27]:

$$\hat{\mathbf{r}}_0 = \arg\max_{\mathbf{r}_0} \frac{\mathbf{G}^H(\mathbf{r}_0)\hat{\mathbf{K}}\mathbf{G}(\mathbf{r}_0)}{|\mathbf{G}(\mathbf{r}_0)|^2} \qquad (11)$$

and

$$\hat{\mathbf{r}}_0 = \arg\max_{\mathbf{r}_0} \frac{1}{\mathbf{G}^H(\mathbf{r}_0)\hat{\mathbf{K}}^{-1}\mathbf{G}(\mathbf{r}_0)} \qquad (12)$$

respectively. It is worth noting that to guarantee the full rank of the sample covariance matrix, the Capon's MFP method needs to use a relatively large number of snapshots (greater than the

rank of the sample covariance matrix) or regularization methods [38]. By contrast, the conventional MFP works even if only one snapshot is available.

### IV. NUMERICAL EXPERIMENTS

In this section, numerical experiments are performed. Section IV-A recalls the principles of SEM that is used in the forward problem to simulate the wave propagation and measurements in a complex environment. Then, sound-source localization using the four-layer model is studied. A summer Mediterranean shallow-water sound-speed profile is considered. The influence of an additional range-dependent variation of the medium (nonflat water bottom and horizontal-sound-speed variation) is quantified in Sections IV-B and IV-C.

#### A. Sound Propagation in a Complex Environment Using SEM

The sound propagation in a complex underwater environment is realized via SEM, which is a time-domain full-wave numerical approach. The accuracy of this numerical method, particularly for sound wave simulation in ocean acoustics, has been validated in the literature [16]. In this section, the wave equation is assumed to be independent of the horizontal azimuth, therefore as mentioned above the 2-D axisymmetric version of the SEM [26] is used. The water depth is $d = 100$ m (in the next section a nonflat bottom will be considered, but its average water depth will be 100 m). A real summer Mediterranean shallow-water sound-speed profile (the data used are from [15, Fig. 1.15]) in the water column is shown by the solid line in Fig. 2; its linear approximation is represented by the dashed line. The sound speed in the bottom is assumed to be $1800 \text{ m} \cdot \text{s}^{-1}$. The density in the water column and ocean bottom is respectively $\rho_1 = 1000 \text{ kg} \cdot \text{m}^{-3}$ and $\rho_2 = 1800 \text{ kg} \cdot \text{m}^{-3}$. The source signal is a Ricker wavelet, i.e., the second derivative of a Gaussian function, with dominant frequency $f_0 = 100$ Hz. To ensure the accuracy of the SEM, the size $\Delta x$ of each element must be smaller than $\Delta x_{\min} = \lambda_{\min} = 6.04$ m, wherein $\lambda_{\min} \approx v_2/(2.5f_0)$ is the minimum typical wavelength for which there is significant energy in the medium. Besides, according to the Courant-Friedrichs-Lewy stability condition of the time-integration scheme, the time step $\Delta t$ should be smaller than $C\Delta x_{\min}/v_3 = 0.0013$ s, where $C$ is the Courant number, which for the SEM is typically 0.4. We thus select an element length $\Delta x = 1$ m (each element includes nine Gauss-Lobatto-Legendre integration points in each direction and thus 81 points in total) and a time step $\Delta t = 2 \times 10^{-5}$ s, which ensure the accuracy of the wave propagation simulation. Moreover, the considered maximum horizontal range is 10 km, thus the total measurement time is set to 5 s such that the direct propagation and the main reflections of sound wave all pass through the measurement stations. By performing a discrete Fourier transform, the sound pressures in the frequency domain can then be obtained.

Fig. 3 compares the sound pressure at the frequency $f = 100$ Hz at $\mathbf{r}'_m = (2000, 0, z'_m)$ m, $z'_m = 5, 10, \dots, 95$, obtained based on the wave number integration method (crosses) and the SEM (circles), in which the sound speed in both cases is the linearized one (i.e., the four-layer model). Both real and
<figure><img src="image_4.png" /><figcaption>Fig. 2. Sound-speed profile in the summer Mediterranean (solid line, the data used are from [15, Fig. 1.15]), and its approximation using a four-layer model (dashed line).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Real (left) and imaginary (right) parts of sound pressure measurements with the four-layer sound-speed profile, obtained from the wave number integration method (crosses) and the SEM (circles). The measurements are performed at a horizontal distance of 2 km from the source, at depths $z'_{m} = 5, 10, \dots, 95$ m.</figcaption></figure>

imaginary parts of the pressure normalized by the corresponding maximum values among the 19 measurements are shown. This example illustrates the accuracy of both the wave number integration method and the SEM. Here, we remark that the computation time of the wave number integration method for obtaining the results in Fig. 3 is very short (especially compared to SEM), approximately a few seconds on a personal

<figure><img src="image_6.png" /><figcaption>Fig. 4. Real (left), imaginary (middle) parts and absolute value (right) of sound pressure measurements obtained from the four-layer sound-speed profile using the wave number integration method (crosses) and the real sound-speed profile using the SEM (circles). The measurements are performed at a horizontal distance of 5 km from the source, at depths $z'_{m} = 5, 10, \dots, 95$ m.</figcaption></figure>

PC. Then, we consider the more realistic sound-speed profile (the solid line in Fig. 2) and the corresponding sound pressures are calculated using the SEM. The circles in Fig. 4 show the normalized real part, imaginary part and norm of the sound pressure at $r'_{m} = (5000, 0, z'_{m})$ m, $z'_{m} = 5, 10, \dots, 95$, as well as the corresponding results (crosses) with sound-speed approximation ($d_{1} = 15$ m, $d_{2} = 70$ m, $c_{1} = 1537.5$ m·s⁻¹, and $c_{2} = 1510$ m·s⁻¹, as shown in Fig. 2) using the wave number integration method. In this case, due to the mismatch of the sound-speed profile between the irregular one and its linearized approximation, the simplified model cannot perfectly replicate the phase of sound pressure but the amplitude is relatively precise. Note that if the data generator for the inverse problem could precisely compute the measured pressure as in Fig. 3, the localization, which essentially matches the measurement to the model by tuning the corresponding model parameters, would be easy. The main purpose of this paper is to evaluate the performance of the localization methods under the model and data uncertainties like the case in Fig. 4. It has been shown in [28] that the MFP approaches are robust with uncertainties since they can efficiently deal with phase error. In the following, the limitations of MFP approaches are tested in more complex and practical scenarios: An additional range-dependent variation of the underwater environment is considered, in which we will study and quantify the source localization accuracy.

## B. Sound-Source Localization With Irregular Bottom Interfaces

In this section, we study the source localization error with respect to bottom interface irregularities. As a reference, the localization error, which is defined by the $L_{2}$ distance between the source location estimate and its actual value

$$e(\hat{r}_{0}) = \|\hat{r}_{0} - r_{0}\|_{2} \qquad (13)$$

with only sound-speed approximation is first plotted. In the forward problem, the data generation is realized via the SEM
<figure><img src="image_7.png" /><figcaption>Fig. 5. Source localization error (absolute error and error relative to wavelength) in shallow water, along with its 95% confidence interval, due to measurement noise and sound-speed linear approximation. SNR varies from -10 dB to 10 dB. The actual sound-speed profile and its approximation follow the solid line and the dash line in Fig. 2.</figcaption></figure>

with the irregular sound-speed profile (the solid line in Fig. 2). Gaussian random noise with various signal-to-noise ratio (SNR) is added to obtain the data for the inverse problem and T = 220 snapshots are used, which guarantees the full rank of the sample covariance matrix. Here, SNR is defined by $20 \log_{10}(\bar{p}/\sigma)$, where $\sigma$ is the standard deviation of the Gaussian random noise and $\bar{p}$ is the average sound pressure. Note that by the Shannon-Nyquist sampling theorem [39], [40], the resolution of leak localization is approximately half of a wavelength. Therefore, throughout this section, both absolute error and error relative to wavelength $\lambda$ is plotted and the error less than a half-wavelength is considered as “acceptable.” The coordinates of measurement stations are $r'_m = (5000, y'_m, z'_m)$ m, $y'_m = -1000, -600, \dots, 1000$ and $z'_m = 5, 15, \dots, 95$. Fig. 5 shows the localization results of conventional and Capon's MFPs with various SNR from -10 dB to 10 dB; the average localization error and the corresponding 95% confidence interval obtained from 30 computations (with different realizations of random noise) are plotted. For high SNRs, the error tends to 0.5 m (around $0.02\lambda$), instead of 0, which can be assumed as the error due to sound-speed approximation. Note that this error is small, which implies that the linear approximation of the sound-speed profile is reasonable in the sense of source localization. As the noise level increases to SNR = -10 dB, the error of both methods increases but Capon's MFP has a relatively higher error. However, both errors are acceptable ($\ll 0.5\lambda$) even for a high noise level SNR = -10 dB.

Then, nonflat water bottom is considered. The data are generated via the SEM with the irregular sound-speed profile and a sinusoidal bottom interface (i.e., range dependent)

$$d(x) = 100 + A_b \sin(\lambda_b x / (2\pi)). \quad (14)$$

Fig. 6 shows the sound field in this case at the time snapshots 0.154, 0.259, 0.609, and 0.756 s, with $A_b = 8$ m and $\lambda_b =$

<figure><img src="image_8.png" /><figcaption>Fig. 6. Snapshots of sound propagation at times 0.154, 0.259, 0.609, and 0.756 s with the summer Mediterranean sound-speed profile (the solid line in Fig. 2) and with a sinusoidal bottom with $\lambda_b = 400$ m and $A_b = 8$ m in (14).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Sound-source localization in the summer Mediterranean using a four-layer model. The water bottom is a sinusoidal function with period $\lambda_b = 1000$ m and amplitude $A_b = 5$ m in (a) and (b) and $A_b = 7$ m in (c) and (d). The source localization methods are conventional MFP in (a) and (c) and Capon MFP in (b) and (d). The distance between the sound source and the measurement plane is 5 km and the frequency is 100 Hz.</figcaption></figure>

400 m in (14). We then perform the source localization assuming the (range independent) four-layer sound-speed profile and flat bottom surface. To evaluate the localization performance of the MFP methods, we first plot the ambiguity surface [28] in the source plane $x = 0$ (by 2-D plotting the output function in (11) or (12) with fixed $x = 0$) in Fig. 7, where $\lambda_b = 1000$ m. When $A_b = 5$ m [Fig. 7(a) and (b)], both methods return an acceptable estimate of sound-source location, although with slight errors. By contrast, when the maximum variation of water depth ($A_b$) is increased to 7 m [Fig. 7(c) and (d)] the source localization is not robust: The reconstructed sound field of conventional MFP exhibits a spurious secondary source and the estimate of Capon's MFP has a large error.

Then, Fig. 8 plots the source localization errors for different amplitudes of the sinusoidal bottom surface, $\lambda_b = 1000$ m in Fig. 8(a), (c), and (e) and $\lambda_b = 400$ m in Fig. 8(b), (d), and (f), and SNR is 20 dB in Fig. 8(a) and (b), 5 dB in Fig. 8(c) and (d), and -5 dB in Fig. 8(e) and (f). For the low noise
<figure><img src="image_10.png" /><figcaption>Fig. 8. Source localization error in shallow water with a sinusoidal bottom. The period of the bottom is 1000 m in (a), (c), and (e) and 400 m in (b), (d), and (f). The SNR is 20 dB in (a) and (b), 5 dB in (c) and (d), and −5 dB in (e) and (f).</figcaption></figure>

level SNR = 20 dB where the influence of noise can be almost neglected, both source localization methods are robust for $A_b \le 6$ m in Fig. 8(a). For the sinusoidal bottom with $\lambda_b = 400$ m, the errors of both methods are larger than in the previous case ($\lambda_b = 1000$ m) due to the presence of a larger number of slopes: The localization error is greater than a half-wavelength when $A_b \ge 2$ m. However, the error of conventional MFP is smaller than the Capon's MFP, which shows that the former method is more robust with respect to modeling error. As SNR decreases (to 5 dB and −5 dB), the localization error increases, except that in Fig. 8(c) and (d) when $A_b = 7$ m, the average error of Capon's MFP decreases. This is because two large local maxima of (12) appear [can be observed as in Fig. 7(d)]: when noise is low, the one corresponding to the spurious secondary source is higher but a higher noise possibly enables the peak corresponding to the actual leak location to be higher. Generally speaking, the influence of noise is relatively small compared to the uncertainty due to nonflat bottom: The threshold of $A_b$ in the sense of error less than a half-wavelength (labeled by the dotted lines in Fig. 8) almost does not change even for a low SNR = −5 dB.

### C. Sound-Source Localization With Horizontally Varying Sound Speed

In this section, a sound speed having a horizontal variation is considered, which may represent some real cases of sound-speed profile due to local variation of temperature and salinity. Note that the sound speed follows the following formula [11]:

<figure><img src="image_11.png" /><figcaption>Fig. 9. Source localization error in shallow water with a range-dependent variation sound speed, using the conventional and Capon's MFP methods. The coefficient of horizontal variation is $h = 250, 500, 750, \text{or } 1000$.</figcaption></figure>

$$c(z, T, S) = 1449.2 + 4.6T - 0.055T^2 + 0.00029T^3 + (1.34 - 0.01T)(S - 35) + 0.016z \quad (15)$$

where $S \in [0, 45]$ is the salinity and $T \in [0, 35^\circ\text{C}]$ is the temperature. In the numerical example in this section, the sound speed is assumed to be $c(x, z) = c_0(z) + (x - 5000)/h$, where $c_0(z)$ is the range-independent sound speed shown by the solid line in Fig. 2. This sound-speed profile $c(x, z)$ can represent a linear variation of salinity or an approximately linear variation of temperature in the horizontal direction. The value $h$ of horizontal variation is set to 250, 500, 750, and then 1000, which corresponds to different levels of horizontal variation of sound speed. The bottom surface is flat and the water depth is $d = 100$ m. Different distances between the sound source and the measurement plane are considered, from 1 to 10 km. Fig. 9 shows the source localization errors using the conventional and Capon's MFP. When the source-receiver plane distance is less than 7 km, a high horizontal variation of sound speed ($h = 250$) results in a large error, since the sound-speed error in this range has a large bias. In the three other cases ($h = 500, 750, \text{and } 1000$), the error is acceptable (almost smaller than a half-wavelength). When the source-receiver plane distance increases to 10 km, i.e., the sound-speed errors in the horizontal direction may compensate each other (the sound speed is overestimated in short range and underestimated in long range), the errors for all the four cases are similar and acceptable. Based on this result, we can conclude that an unbiasedly horizontal-sound-speed error (the average error in the horizontal range is 0) does not significantly affect the sound localization result. Furthermore, the error of Capon's MFP is slightly higher than the conventional MFP, which is similar to the previous result.

### V. CONCLUSION

In shallow-water sound-source localization problems, the sound-speed profile is always simplified such that the sound propagation can be rapidly computed. In this paper, a four-layer
model was considered in which the medium was assumed to be horizontally stratified and vertically divided into four layers, and in each layer the sound speed was constant or linear. Based on this assumption, the wave propagation can be easily computed via a onefold numerical integration. This paper has evaluated the suitability of this model in terms of source localization in a realistic environment. The acoustical measurements were simulated via a time-domain full wave simulation technique called the SEM. The vertical sound speed followed a summer Mediterranean shallow-water profile. Two cases of range-dependent variation of the underwater environment were also considered. First, when the bottom was assumed to be non-flat, the source localization accuracy with respect to different levels of bottom irregularity was evaluated. The threshold for bottom irregularity parameters for accurate source localization (error less than a half-wavelength) with various noise level was shown. Second, a horizontally varying sound speed that may represent variation of temperature and salinity was considered, which showed that the model simplification of sound-speed profile is reasonable if it is approximately unbiased in the horizontal direction.

The numerical results in this paper show two types of range-dependent underwater environments: seabed irregularity (slopes) and horizontal variation in the water column due to temperature, salinity, etc. However, real shallow-water environments are much more complex than these examples and should be studied in the future. For example, very different seafloors (gravels, boulders, or very fine mud) can be found within a propagation range of a few kilometers. Besides, the seabed may be more reasonably assumed to be made of many layers with nonparallel interfaces. Generally speaking, to simulate the wave propagation in these environments, the axisymmetric 2-D SEM cannot be used but the 3-D version of SEM is necessary, which requires a higher computation power.

# APPENDIX A
## SOLUTION OF THE DEPTH-SEPARATED WAVE EQUATION

The depth-separated wave equation (8) with the four-layer model has a general solution [15], [19]:

$$ \psi(k_r, z) = \psi_g(k_r, z) + \psi_s(k_r, z) \quad (A.1) $$

where $\psi_g(k_r, z)$ is the solution of the homogeneous depth-separated wave equation [i.e., when the right-hand side of (8) is 0]

$$ \psi_g(k_r, z) = \begin{cases} A_1^+(k_r)e^{ik_{z,1}z} + A_1^-(k_r)e^{-ik_{z,1}z}, & z \in [0, d_1) \\ A_2^+(k_r)\mathrm{Ai}(\xi) + A_2^-(k_r)[\mathrm{Ai}(\xi) - \mathrm{Bi}(\xi)], & z \in [d_1, d_2) \\ A_3^+(k_r)e^{ik_{z,2}z} + A_3^-(k_r)e^{-ik_{z,2}z}, & z \in [d_2, d) \\ A_4^+(k_r)e^{ik_{z,3}(z-d)}, & z \in [d, \infty) \end{cases} \quad (A.2) $$

and

$$ \begin{aligned} \psi_s(k_r, z) = & \begin{cases} S_f \frac{e^{ik_{z,1}|z-z_0|}}{4\pi i k_{z,1}} \mathbf{1}_{[0,d_1)}(z), & \text{if } z_s \in [0, d_1) \\ \frac{S_f}{4\pi} f_-(k_r, z) \mathbf{1}_{[d_1, z_s)}(z) + \frac{S_f}{4\pi} f_+(k_r, z) \mathbf{1}_{[z_s, d_2)}(z), & \text{if } z_s \in [d_1, d_2) \\ S_f \frac{e^{ik_{z,2}|z-z_0|}}{4\pi i k_{z,2}} \mathbf{1}_{[d_2, d)}(z), & \text{if } z_s \in [d_2, d). \end{cases} \end{aligned} \quad (A.3) $$

Here, 1 is the indicator function, the vertical wave numbers are

$$ k_{z,i} = \begin{cases} \sqrt{k_i^2 - k_r^2}, & |k_r| \le k_i \\ i\sqrt{k_r^2 - k_i^2}, & |k_r| > k_i \end{cases}, \quad k_i = \frac{2\pi f}{c_i}, \quad i = 1, 2, 3 \quad (A.4) $$

$$ f_-(k_r, z) = \frac{2(\omega^2 a)^{-1/3} [\mathrm{Ai}(\xi_0) - \mathrm{Bi}(\xi_0)] \mathrm{Ai}(\xi)}{\mathrm{Ai}'(\xi_0)[\mathrm{Ai}(\xi_0) - \mathrm{Bi}(\xi_0)] - \mathrm{Ai}(\xi_0)[\mathrm{Ai}'(\xi_0) - \mathrm{Bi}'(\xi_0)]} \quad (A.5) $$

$$ f_+(k_r, z) = \frac{2(\omega^2 a)^{-1/3} \mathrm{Ai}(\xi_0) [\mathrm{Ai}(\xi) - \mathrm{Bi}(\xi)]}{\mathrm{Ai}'(\xi_0)[\mathrm{Ai}(\xi_0) - \mathrm{Bi}(\xi_0)] - \mathrm{Ai}(\xi_0)[\mathrm{Ai}'(\xi_0) - \mathrm{Bi}'(\xi_0)]} \quad (A.6) $$

where $\mathrm{Ai}$ and $\mathrm{Bi}$ are Airy functions of the first and second kinds, $\xi$ is the variable transformation

$$ \xi = (\omega^2 a)^{-2/3} [k_r^2 - \omega^2 (az + b)] \quad (A.7) $$

and $\xi_0$ is obtained by inserting $z_0$ into (A.7). The coefficients $A_1^+, A_1^-, A_2^+, A_2^-, A_3^+, A_3^-$, and $A_4^+$ can be computed by enforcing the pressure-free boundary at $z=0$, i.e.,

$$ \psi_1(k_r, 0) = 0 \quad (A.8) $$

and the continuity of vertical displacement and pressure at $z = d_1, d_2, d$, i.e.,

$$ \frac{\partial \psi_1}{\partial z}(k_r, \delta) = \frac{\partial \psi_2}{\partial z}(k_r, \delta) \quad (A.9) $$

and

$$ \psi_1(k_r, \delta) = \psi_2(k_r, \delta) \quad (A.10) $$

where $\delta = d_1, d_2, d$. The solution can be obtained by solving the linear equation

$$ S_i = QA \quad (A.11) $$

if the source is in the $i$th layer ($i=1, 2, 3$). Then, in that equation, we have (A.12) and (A.13), shown at the top of the next page, where $C_i = A_i - iB_i$, and then
$$ \mathbf{A} = (A_1^+, A_1^-, A_2^+, A_2^-, A_3^+, A_3^-, A_4^+)^\top \quad (A.12) $$

$$ \mathbf{Q} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ i k_{z,1} e^{i k_{z,1} d_1} & -i k_{z,1} e^{-i k_{z,1} d_1} & (\omega^2 a)^{\frac{1}{3}} \mathrm{Ai}'(\xi_{d_1}) & (\omega^2 a)^{\frac{1}{3}} \mathrm{Ci}'(\xi_{d_1}) & 0 & 0 & 0 & 0 \\ e^{i k_{z,1} d_1} & e^{-i k_{z,1} d_1} & -\mathrm{Ai}(\xi_{d_1}) & -\mathrm{Ci}(\xi_{d_1}) & 0 & 0 & 0 & 0 \\ 0 & 0 & (\omega^2 a)^{\frac{1}{3}} \mathrm{Ai}'(\xi_{d_2}) & (\omega^2 a)^{\frac{1}{3}} \mathrm{Ci}'(\xi_{d_2}) & i k_{z,2} e^{i k_{z,2} d_2} & -i k_{z,2} e^{-i k_{z,2} d_2} & 0 & 0 \\ 0 & 0 & \mathrm{Ai}(\xi_{d_2}) & \mathrm{Ci}(\xi_{d_2}) & -e^{i k_{z,2} d_2} & -e^{-i k_{z,2} d_2} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & i k_{z,2} e^{i k_{z,2} d} & -i k_{z,2} e^{-i k_{z,2} d} & -i k_{z,3} \\ 0 & 0 & 0 & 0 & 0 & \rho_1 e^{i k_{z,2} d} & \rho_1 e^{-i k_{z,2} d} & -\rho_2 \end{pmatrix} \quad (A.13) $$

$$ \mathbf{S}_1 = -\frac{S_f}{4\pi} \left( 0, (\omega^2 a)^{\frac{1}{3}} \frac{\partial f_-}{\partial \xi}(k_r, d_1), -f_-(k_r, d_1), (\omega^2 a)^{\frac{1}{3}} \frac{\partial f_+}{\partial \xi}(k_r, d_2), f_+(k_r, d_2), 0, 0 \right)^\top \quad (A.14) $$

$$ \mathbf{S}_2 = -\frac{S_f}{4\pi} \left( \frac{e^{ik_{z,1}z_0}}{ik_{z,1}}, e^{ik_{z,1}(d_1-z_0)}, \frac{e^{ik_{z,1}(d_1-z_0)}}{ik_{z,1}}, 0, 0, 0, 0 \right)^\top \quad (A.15) $$

and

$$ \mathbf{S}_3 = \frac{S_f}{4\pi} \left( 0, 0, 0, e^{ik_{z,2}(z_0-d_2)}, \frac{e^{ik_{z,2}(z_0-d_2)}}{ik_{z,2}}, -e^{ik_{z,2}(d-z_0)}, -\rho_1 \frac{e^{ik_{z,2}(d-z_0)}}{ik_{z,2}} \right)^\top . \quad (A.16) $$

## APPENDIX B
### WAVENUMBER INTEGRATION

Let us briefly recall how the Helmholtz equation (5) can be obtained by numerically computing the inverse Hankel transform [15], [19], [20]. First, the integration domain in (9) is truncated. It can be observed that the integral kernel decays rapidly after the maximum wave number in the medium $k_2 = 2\pi f/c_2$ and quickly reaches a certain point, denoted as $k_{max}$, at which the matrix $\mathbf{Q}$ is singular. The numerical integral is thus truncated right before $k_{max}$.

The depth-separated wave equation $\psi(k_r, z)$ contains discrete singularities in the interval $[k_3, k_2]$. For this reason, a complex contour integration method [15], based on the Cauchy's residue theorem, is used to obtain a numerically stable integration: The integral domain is moved from the real axis (from 0 to $k_{max}$) to the complex plane with a small vertical displacement $\epsilon$ ($\epsilon \ll k_{max}$):

$$ \psi(r, z) \approx \int_{0}^{k_{max}} \psi((k_r - i\epsilon), z) J_0((k_r - i\epsilon)r)(k_r - i\epsilon)dk_r. \quad (B.1) $$

The right-hand side of (B.1) guarantees the nonsingularity of the integration and smooths the integral kernel.

The integral kernel of (B.1) oscillates more significantly in the interval $[k_3, k_2]$ than in the other two regions $[0, k_3)$ and

$(k_2, k_{max})$. Therefore, a smaller discretization spacing in the interval $[k_3, k_2]$ (denoted by $\Delta_2$) than the other two (the spacing is $\Delta_1$ and $\Delta_3$) is well advised. Thus, the Helmholtz equation at each point with range $r$ and depth $z$ can be computed via the numerical integration of (B.1):

$$ \hat{\psi}(r, z) = \left( \Delta_1 \sum_{l=1}^{L_1} + \Delta_2 \sum_{l=L_1+1}^{L_1+L_2} + \Delta_3 \sum_{l=L_1+L_2+1}^{L_1+L_2+L_3} \right) \times \psi((k_l - i\epsilon), z) J_0((k_l - i\epsilon)r)(k_l - i\epsilon) \quad (B.2) $$

where $k_l$ s are the discrete points of $k_r$, $L_1 = k_3/\Delta_1$, $L_2 = (k_2 - k_3)/\Delta_2$, and $L_3 = (k_{max} - k_2)/\Delta_3$.

## ACKNOWLEDGMENT

The authors would like to thank P. Cristini for discussion and two anonymous reviewers for their valuable comments and suggestions, which greatly improved this paper.

## REFERENCES

[1] W. A. Kuperman and J. F. Lynch, "Shallow-water acoustics," Phys. Today, pp. 55-61, Oct. 2004.

[2] J. Tabrikian and J. L. Krolik, "Robust maximum-likelihood source localization in an uncertain shallow-water waveguide," J. Acoust. Soc. Amer., vol. 101, no. 1, pp. 241-249, 1997.

[3] J. Tabrikian and J. L. Krolik, "Barankin bounds for source localization in an uncertain ocean environment," IEEE Trans. Signal Process., vol. 47, no. 11, pp. 2917-2927, Nov. 1999.

[4] R. M. Hamson and R. M. Heitmeyer, "Environmental and system effects on source localization in shallow water by the matched-field processing of a vertical array," J. Acoust. Soc. Amer., vol. 86, no. 5, pp. 1950-1959, 1989.

[5] S. E. Dosso, "Environmental uncertainty in ocean acoustic source localization," Inverse Probl., vol. 19, pp. 419-431, 2003.

[6] S. Finette, "A stochastic representation of environmental uncertainty and its coupling to acoustic wave propagation," J. Acoust. Soc. Amer., vol. 120, no. 5, pp. 2567-2579, 2006.

[7] X. Wang, S. Khazaei, and P. Sagaut, "Sound source localization in a randomly inhomogeneous medium using matched statistical moment method," J. Acoust. Soc. Amer., vol. 138, no. 6, pp. 3896-3906, 2015.

[8] W. Xu and H. Schmidt, "System-orthogonal functions for sound speed profile perturbation," IEEE J. Ocean. Eng., vol. 31, no. 1, pp. 156-169, Jan. 2006.

[9] M. D. Collins and E. K. Westwood, "A higher-order energy-conserving parabolic equation for range-dependent ocean depth, sound speed, and density," J. Acoust. Soc. Amer., vol. 89, no. 3, pp. 1068-1075, 1991.

[10] F. B. Jensen and C. M. Ferla, "Numerical solutions of range-dependent benchmark problems in ocean acoustics," J. Acoust. Soc. Amer., vol. 87, no. 4, pp. 1499-1510, 1990.

[11] H. Medwin, "Speed of sound in water: A simple equation for realistic parameters," J. Acoust. Soc. Amer., vol. 58, no. 6, pp. 1318-1319, 1975.
[12] N. R. Chapman, S. Chin-Bing, D. King, and R. B. Evans, "Benchmarking geoo acoustic inversion methods for range-dependent waveguides," *IEEE J. Ocean. Eng.*, vol. 28, no. 3, pp. 320–330, Jul. 2003.

[13] P. Abbot, S. Celuzza, I. Dyer, B. Gomes, J. Fulford, and J. Lynch, "Effects of east china sea shallow-water environment on acoustic propagation," *IEEE J. Ocean. Eng.*, vol. 28, no. 2, pp. 192–211, Apr. 2003.

[14] G. F. Edelmann, T. Akal, W. S. Hodgkiss, S. Kim, W. A. Kuperman, and H. C. Song, "An initial demonstration of underwater acoustic communica- tion using time reversal," *IEEE J. Ocean. Eng.*, vol. 27, no. 3, pp. 602–609, Jul. 2002.

[15] F. B. Jensen, W. A. Kuperman, M. B. Porter, and H. Schmidt, *Computational Ocean Acoustics*, 2nd ed. Berlin, Germany: Springer-Verlag, 2011.

[16] P. Cristini and D. Komatitsch, "Some illustrative examples of the use of a spectral-element method in ocean acoustics," *J. Acoust. Soc. Amer.*, vol. 131, no. 3, pp. EL229-EL235, 2012.

[17] R. B. Evans, "A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom," *J. Acoust. Soc. Amer.*, vol. 74, no. 1, pp. 188-195, 1983.

[18] M. D. Collins, "A split-step Padé solution for the parabolic equation method," *J. Acoust. Soc. Amer.*, vol. 93, no. 4, pp. 1736-1742, 1993.

[19] F. R. DiNapoli and R. L. Deavenport, "Theoretical and numerical Green's function field solution in a plane multilayered medium," *J. Acoust. Soc. Amer.*, vol. 67, no. 1, pp. 92-105, 1980.

[20] H. Schmidt and F. B. Jensen, "A full wave solution for propagation in mul-tilayered viscoelastic media with application to Gaussian beam reflection at fluid-solid interfaces," *J. Acoust. Soc. Amer.*, vol. 77, no. 3, pp. 813-825, 1985.

[21] E. K. Westwood, C. T. Tindle, and N. R. Chapman, "A normal mode model for acousto-elastic ocean environments," *J. Acoust. Soc. Amer.*, vol. 100, no. 6, pp. 3631-3645, 1996.

[22] S. J. Levinson, E. K. Westwood, R. A. Koch, S. K. Mitchell, and C. V. Sheppard, "An efficient and robust method for underwater acous- tic normal-mode computations," *J. Acoust. Soc. Amer.*, vol. 97, no. 3, pp. 1576-1585, 1995.

[23] C. A. Clark and K. B. Smith, "An efficient normal mode solution to wave propagation prediction," *IEEE J. Ocean. Eng.*, vol. 33, no. 4, pp. 462-476, 2008.

[24] A. T. Patera, "A spectral element method for fluid dynamics: Laminar flow in a channel expansion," *J. Comput. Phys.*, vol. 54, no. 3, pp. 468-488, 1984.

[25] D. Komatitsch and J.-P. Vilotte, "The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures," *Bull. Seismol. Soc. Amer.*, vol. 88, no. 2, pp. 368-392, 1998.

[26] A. Bottero, P. Cristini, D. Komatitsch, and M. Asch, "An axisymmet- ric time-domain spectral-element method for full-wave simulations: Ap- plication to ocean acoustics," *J. Acoust. Soc. Amer.*, vol. 140, no. 5, pp. 3520-3530, 2016.

[27] A. B. Baggeroer, W. A. Kuperman, and P. N. Mikhalevsky, "An overview of matched field methods in ocean acoustics," *IEEE J. Ocean. Eng.*, vol. 18, no. 4, pp. 401-424, Oct. 1993.

[28] A. Tolstoy, *Matched Field Processing for Underwater Acoustics*. Singapore: World Scientific, 1993.

[29] M. D. Collins and W. A. Kuperman, "Inverse problems in ocean acoustics," *Inverse Probl.*, vol. 10, no. 5, pp. 1023-1040, 1994.

[30] Z.-H. Michalopoulou and M. B. Porter, "Matched-field processing for broad-band source localization," *IEEE J. Ocean. Eng.*, vol. 21, no. 4, pp. 384-392, Oct. 1996.

[31] H. P. Bucker, "Use of calculated sound fields and matched field detection to locate sound sources in shallow water," *J. Acoust. Soc. Amer.*, vol. 59, no. 2, pp. 368-373, 1976.

[32] R. G. Fizell and S. C. Wales, "Source localization in range and depth in an Arctic environment," *J. Acoust. Soc. Amer.*, vol. 78, no. S1, pp. S57-S58, 1985.

[33] H. Krim and M. Viberg, "Two decades of array signal processing research," *IEEE Signal Process. Mag.*, vol. 13, no. 4, pp. 67-94, Jul. 1996.

[34] X. Wang, B. Quost, J.-D. Chazot, and J. Antoni, "Estimation of multiple sound sources with data and model uncertainties using the EM and eviden- tial EM algorithms," *Mech. Syst. Signal Process.*, vol. 66-67, pp. 159-177, 2016.

[35] X. Wang, B. Quost, J.-D. Chazot, and J. Antoni, "Iterative beamforming for identification of multiple broadband sound sources," *J. Sound Vib.*, vol. 365, pp. 260-275, 2016.

[36] J. Capon, "High-resolution frequency-wavenumber spectrum analysis," *Proc. IEEE*, vol. 57, no. 8, pp. 1408-1418, Aug. 1969.

[37] S. D. Somasundaram and N. H. Parsons, "Evaluation of robust Capon beamforming for passive sonar," *IEEE J. Ocean. Eng.*, vol. 36, no. 4, pp. 686-695, Oct. 2011.

[38] O. Ledoit and M. Wolf, "A well-conditioned estimator for large-dimensional covariance matrices," *J. Multivariate Anal.*, vol. 88, no. 2, pp. 365-411, 2004.

[39] C. E. Shannon, "Communication in the presence of noise," *Proc. IRE*, vol. 37, no. 1, pp. 10-21, 1949.

[40] A. J. Jerri, "The Shannon sampling theorem—its various extensions and applications: A tutorial review," *Proc. IEEE*, vol. 65, no. 11, pp. 1565-1596, Nov. 1977.