# Sound source localization in a randomly inhomogeneous medium using matched statistical moment method

Xun Wang, Shahram Khazaie, Pierre Sagaut

► To cite this version:

Xun Wang, Shahram Khazaie, Pierre Sagaut. Sound source localization in a randomly inhomogeneous medium using matched statistical moment method. Journal of the Acoustical Society of America, 2015, 138 (6), pp.3896. 10.1121/1.4938238. hal-01276517

HAL Id: hal-01276517

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

Submitted on 26 Feb 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Sound source localization in a randomly inhomogeneous medium using matched statistical moment method

Xun Wang,ª) Shahram Khazaie, and Pierre Sagaut

Aix Marseille Université, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451 Marseille, France

This paper investigates the problem of sound source localization from acoustical measurements obtained by an array of microphones. The sound propagation medium is assumed to be randomly inhomogeneous, being modelled by a random function of space. In this case, classical source local-ization methods (e.g., beamforming, near-field acoustical holography, and time reversal) cannot be used anymore. Therefore, an approach based on the statistical moments of acoustical measurement is proposed to solve the aforementioned problem. In this work, a Karhunen Loève expansion is used so that the random medium can be represented by a small number of uncorrelated and identi-cally distributed random variables. The statistical characteristics of the measurements in terms of probability density function and statistical moments are also studied. Then, the sound source is localized by minimizing the error of statistical moments between the real measurements obtained from the microphone array and the measurements simulated from an assumed source. Finally, a numerical example is introduced to justify the proposed method. This experiment shows that the random field can be replicated by a very small number of random variables, the statistical moments of measurements guarantee the convergence, and the source location can be accurately estimated using the proposed source localization method.

## I. INTRODUCTION

Source localization is an important issue in acoustics engineering, seismology, and oceanic engineering. Beamforming,¹,² near-field acoustical holography (NAH),³⁻⁵ and time reversal⁶⁻¹¹ are all classical methods dedicated to the aforementioned issue. Beamforming estimates the arrival direction of plane wave or the location of point source via the delay of signal arrival measured by the microphones. NAH back-propagates the sound field in the wavenumber domain from the measurement plane to the source plane, which ensures a high spatial resolution by taking evanescent waves into account. Time reversal method reverses the received signal in time and retransmits it back into the medium, which can finally refocus on the source.

However, contrary to the ideal environment assumed in the above approaches, uncertainties may pervade the sound propagation and measurement process. In particular, the uncertainty of sound propagation medium may significantly affect the source localization, for instance in the ocean environment due to limited knowledge of physical properties of the medium,¹² or in an environment with successive and irregular flows that bring in strong randomness.¹³ The effect of medium uncertainties in the sound propagation and source localization has been discussed in the literature.¹²⁻¹⁸ In these works, the uncertainty of medium is described by one or a few parameters, e.g., a random but uniform sound speed in the space. However, in some real cases, the situation may be more complicated: The physical properties at every point in

the medium may be different and random. In this paper, the sound source localization problem is considered under the assumption of a randomly inhomogeneous medium, in which the acoustical field cannot be precisely measured, but only its statistical information is available. More specifically, the sound propagation medium is modelled by a random field, whose expectation and covariance function are predefined. In this case, the classical approaches cannot be used anymore. Beamforming and NAH are based on the assumption of deterministic and homogeneous medium. Time reversal works better in an inhomogeneous medium with multiple scatterings, reflections, and refraction than in the homogeneous environment in the sense of super-resolution.⁶ However, it still requires one to know the precise sound speed at every point in the medium in the back-propagation process. Therefore, this paper proposes a new source localization approach to deal with the case of a randomly inhomogeneous medium.

In order to investigate the uncertainties of sound propagation and measurements, Monte Carlo simulation can be used by repeatedly sampling from the random sound field. However, the complexity of the random medium (large number of dependent random variables) prevents one from solving the problem. In order to overcome this difficulty, the representation of the random field has to be simplified. The spectral representation method¹⁹ is able to represent a Gaussian random field as a summation of trigonometric functions with random amplitudes and phases. This approach is not optimal as far as computational costs are concerned. By using fast Fourier transform, the latter can be dramatically decreased.²⁰ However, this approach is still limited to

<footnote>ª) Electronic mail: xun.wang@univ amu.fr</footnote>
the assumption of Gaussian distribution so that non-Gaussian random field cannot be represented. Polynomial chaos expansion<sup>21–23</sup> is a way of representing an output random variable with a given distribution (not restricted to Gaussian) as a function of a few random variables in the form of polynomial expansion, but it cannot be used to represent the (input) random field in the case of this paper. Karhunen Loève (K L) expansion<sup>21,24,25</sup> can efficiently represent a random field via a small number of random variables. The considered random field is not restricted to Gaussian, but for any second-order field. In this approach, the expectation and covariance function of the random field are the only information needed to be predefined. This assumption is reasonable for many cases including the random sound field considered in this paper.

Given a realization of random medium, the sound propagation and measurement can be simulated. In some special cases of deterministic medium, the Green’s function of wave equation can be analytically computed,<sup>12,26–28</sup> such that the sound field at every point of the medium can be easily obtained. However, for a general inhomogeneous medium (the case of this paper), some numerical methods have to be employed, such as finite difference, finite elements, boundary elements, spectral elements, etc. The dimension of the medium, the complexity of its boundaries, the fluctuation level of the heterogeneities, the Sommerfeld radiation condition at infinity (for unbounded media), and the linear (small deformations) or nonlinear (large deformations) behavior of the underlying medium are some of the influential factors on the efficiency and accuracy of these numerical methods. A detailed review and comparison between these methods is out of the scope of this presentation and can be found, for instance, in Ref. 29. Spectral element method (SEM) is a high-order variational discretization for partial differential equations based on the original ideas of Ref. 30 for the Chebyshev polynomials. It is then developed by Ref. 31 for the Legendre polynomial basis. This method results in a diagonal representation of the mass matrix which significantly reduces the computational cost. As in the higher order finite element method, numerical wave dispersion, which depends highly on the order of finite element shape function, is reduced. Applications of the SEM in two dimensional (2D) and three dimensional (3D) elastodynamics have shown high accuracies and weak numerical dispersions.<sup>32–34</sup> Hence, this method is increasingly used in multidimensional problems in acoustics and elastodynamics notably due to its exponential convergence rate.<sup>35</sup> The SEM can be related to the hp-finite element method,<sup>36</sup> in which the mesh is refined using a suitable combination of h- and p-refinements (subdivision of the elements into smaller ones and increasing their polynomial degree) that result in an exponentially convergent method. The numerical code used in this study is developed based on the SEM. By combining the K L expansion and SEM, propagation of the input uncertainty of random medium to the output uncertainty of acoustical measurement is achieved. Then, statistical properties, e.g., probability density function (PDF) and statistical moments (expectation and covariance matrix), of the acoustical field can be obtained.

In this work, the expectation and covariance matrix of acoustical measurements are the data for estimating the sound source location. Each measurement is obtained from a microphone array and assumed to be a realization of the random output of sound pressure. The estimated source location is decided by a minimum error criterion of statistical moment: the point minimizing the error of statistical moments between its simulated measurement and the real measurement is considered as the source estimate. The organization of this paper is as follows. Section II begins with a description of the model. Then, the K L expansion is used to simplify the representation of the random field of sound propagation medium. The statistical properties of acoustical measurement are also investigated. The matched statistical moment method for sound source localization is introduced in Sec. III; the estimation equations are explicitly given. Section IV presents experimental results on simulated data. In this experiment, the K L expansion can efficiently represent the random field via a small number of random variables. Furthermore, the statistical moments of measurement can converge within a reasonable sample size. Finally, the source localization methods using both expectation and covariance matrix errors are proved to be effective.

## II. STOCHASTIC RESPONSE SURFACE OF SOUND PROPAGATION IN A RANDOMLY INHOMOGENEOUS MEDIUM

In this section, the sound propagation model is introduced first. The K-L expansion is then used to represent the random field. Finally, the statistical characteristics of acoustical measurement are investigated.

### A. Model description

In this work, a randomly inhomogeneous sound propagation medium in a region $\Gamma$ is considered. Let $\{\xi_r : r \in \Gamma\}$ be a second-order random field defined on $\Gamma$ with covariance function $C(r_1, r_2)$, where $r_1$ and $r_2$ are any two points in $\Gamma$. The sound speed at each point $r$ is a function of the random variable $\xi_r$, denoted by $v(r, \xi_r)$. It is assumed that the sound field is produced by a sound source located at $r_s$. The sound propagation process follows the wave equation; the sound pressure $\tilde{p}(r, t)$ at position $r$ and time $t$ is governed by

$$ \nabla^2 - \frac{1}{v^2(r, \xi_r)} \frac{\partial^2}{\partial t^2} \tilde{p}(r, t) = -\delta(r - r_s) \tilde{W}(t), \quad (1) $$

in which $\tilde{W}(t)$ is a function of source signal in the time domain and $\delta$ is the Dirac delta function. Applying Fourier transform to the both sides of Eq. (1) with respect to $t$ results in the following Helmholtz equation for the sound field $p(r, f)$ in the frequency domain:

$$ [\nabla^2 + k^2(r, \xi_r)]p(r, f) = -\delta(r - r_s)W(f), \quad (2) $$

in which $k(r, \xi_r) = (2\pi f)/[v(r, \xi_r)]$ is the wavenumber, $f$ is the frequency, and $W(f)$ is the Fourier transform of $\tilde{W}(t)$. Furthermore, $M$ microphones are placed at discrete locations $r'_m, m = 1, \dots, M$. The main purpose of this paper is to use
the sound pressures measured by the microphone array to localize the sound source. Figure 1 shows a 2D example of the proposed model. The colours in this figure indicate a realization of a randomly inhomogeneous sound speed field, the black point and crosses stand for the sound source and microphones, respectively.

## B. Random medium representation

Since the random field is defined at each point $\mathbf{r} \in \Gamma$, an infinite (and uncountable) number of random variables is needed for its complete description. In this section, the K L expansion is proposed to simplify the random field representation via a small number of random variables. The K L expansion of the random field $\xi_{\mathbf{r}}$ reads

$$ \xi_{\mathbf{r}} = \xi_{\mathbf{r}} + \sum_{n=1}^{\infty} \sqrt{\lambda_n \eta_n f_n(\mathbf{r})}, \quad (3) $$

where $\xi_{\mathbf{r}} = \mathbb{E}(\xi_{\mathbf{r}})$ is the expectation function of the random field, $\lambda_n (\lambda_1 \ge \lambda_2 \ge \dots)$ and $f_n(\mathbf{r})$ are the eigenvalues and eigenfunctions of the covariance function $C(\mathbf{r}_1, \mathbf{r}_2)$, and $\eta_n (n \in \mathbb{N})$ are a series of uncorrelated random variables. Here the eigenvalues and eigenfunctions are defined as the solutions of the Fredholm's integral equation of second kind

$$ \int_{\Gamma} C(\mathbf{r}_1, \mathbf{r}_2) f_n(\mathbf{r}_2) d\mathbf{r}_2 = \lambda_n f_n(\mathbf{r}_1). \quad (4) $$

If $\{\xi_{\mathbf{r}} : \mathbf{r} \in \Gamma\}$ is assumed as a Gaussian random field, $\eta_n (n \in \mathbb{N})$ are independent standard Gaussian random variables. If the covariance function is further assumed as $C(\mathbf{r}_1, \mathbf{r}_2) = \exp[-(\|\mathbf{r}_1 - \mathbf{r}_2\|)/\ell_c]$, $\lambda_n$ and $f_n(\mathbf{r})$ in Eq. (4) can be analytically solved (cf. Sec. 2.1 of Ref. 37). Here, $\ell_c$ is the correlation length of the random field. Moreover, the series in Eq. (3) can be approximated by a finite number of terms

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Sound propagation and measurement in a randomly inhomogeneous medium. The point and crosses stand for the sound source and microphones, respectively.</figcaption></figure>

$$ \hat{\xi}_{\mathbf{r}} = \xi_{\mathbf{r}} + \sum_{n=1}^{N} \sqrt{\lambda_n \eta_n f_n(\mathbf{r})}. \quad (5) $$

Note that in this approximation, those terms with small eigenvalues (playing a relatively weak role in the random field representation) are not taken into account. The global mean square error of the approximation can be computed as$^{24}$

$$ \int_{\Gamma} \mathbb{E}(\hat{\xi}_{\mathbf{r}} - \xi_{\mathbf{r}})^2 d\mathbf{r} = \int_{\Gamma} \sigma^2(\mathbf{r}) d\mathbf{r} - \sum_{n=1}^{N} \lambda_n, \quad (6) $$

in which $\sigma^2(\mathbf{r})$ is the local variance of the random field. The number $N$ of terms in Eq. (5) can be decided by letting the error [Eq. (6)] less than a given threshold.

In the numerical calculation, the considered region $\Gamma$ is discretized by a grid with a finite number of points, denoted as $\Gamma^\# = (\mathbf{r}_1, \dots, \mathbf{r}_I)^\mathrm{T}$. Let $\Xi = (\xi_{\mathbf{r}_1}, \dots, \xi_{\mathbf{r}_I})^\mathrm{T}$, $\lambda_n (\lambda_1 \ge \lambda_2 \ge \dots)$, and $f_n(\Gamma^\#) (n = 1, \dots, I)$ be the eigenvalues and eigenvectors of the covariance matrix of $\Xi$. Then, the $N$-truncated (keeping the most significant modes in terms of large eigenvalues) K L expansion of the random field in the grid is

$$ \hat{\Xi} = \mathbb{E}(\Xi) + \sum_{n=1}^{N} \sqrt{\lambda_n \eta_n f_n(\Gamma^\#)}. \quad (7) $$

It should be pointed out that Eq. (7) enables the random sound field to be easily replicated using a small number of uncorrelated and identically distributed random variables.

## C. Statistical moments of measurement

In this section, the statistical characteristics in the form of PDF and statistical moments are computed. By sampling $R$ realizations from the random variables $\eta_1, \dots, \eta_N$ and substituting them back into Eq. (7) [or Eq. (5) in the continuous case], $R$ realizations of the random sound speed field $v(\mathbf{r}_i, \xi_{\mathbf{r}_i})$, $i = 1, \dots, I$ (or in the continuous case $v(\mathbf{r}, \xi_{\mathbf{r}})$, $\mathbf{r} \in \Gamma$), can be obtained.

Note that the sound field is inhomogeneous, in which case the wave Eqs. (1) and (2) have no analytical solution. Therefore, the SEM is used to compute the sound propagation process; the sound pressure at any point in the grid can be computed. For the sake of source localization, $R$ realizations of sound pressures in the frequency domain at the positions of microphones $\mathbf{r}'_m$, $m = 1, \dots, M$, are computed. Then, the estimated PDFs, as well as the sample statistical moments (for example, the expectation and covariance matrix), of the measurements can be computed. In Sec. III, the statistical moments will be used to localize the sound source.

## III. SOUND SOURCE LOCALIZATION USING MATCHED STATISTICAL MOMENT METHOD

The aim of this section is to use the sound pressures measured by the microphone array to localize the unknown sound source. At a given frequency $f$, the $r$th measurement from the $m$th microphone is denoted by $p_{mr}$ ($m = 1, \dots, M$, $r = 1, \dots, R$). Then, the statistical moments of the measurement can be computed: the sample first-order raw
(expectation) and second-order central (covariance matrix) moments are obtained by

$$M_1 = \frac{1}{R} \sum_{r=1}^{R} (|p_{1r}|, \dots, |p_{Mr}|)^\mathrm{T} \quad (8)$$

and

$$M_2 = \frac{1}{R-1} \sum_{r=1}^{R} (p_{1r} - p_1, \dots, p_{Mr} - p_M) \times (p_{1r} - p_1, \dots, p_{Mr} - p_M)^\mathrm{H}, \quad (9)$$

in which $(\cdot)^\mathrm{H}$ is the operation of vector conjugate transpose and $p_m = (1/R) \sum_{r=1}^{R} p_{mr}$.

Let $\mathbf{r}_s$ denote the real position of sound source. Furthermore, prior knowledge that the real position belongs to a region $\Gamma_1$ is available, i.e., $\mathbf{r}_s \in \Gamma_1 \subset \Gamma$. In most cases, such knowledge can be offered. For instance, the plane of the source is known or predefined, which is the necessary information for using classical beamforming and NAH. Then, the region $\Gamma_1$ is discretized by a grid $\Gamma_1^\#$ with $H$ points. For each $\mathbf{r}_h \in \Gamma_1^\# (h = 1, \dots, H)$, $R'$ measurements with different realizations of random medium are simulated using the method proposed in Sec. II denoted as $\mathbf{p}(h) = \{p_{mr}(h) : m = 1, \dots, M, r = 1, \dots, R'\}$. Similarly, the corresponding first- and second-order moments are obtained by

$$M_1(h) = \frac{1}{R'} \sum_{r=1}^{R'} (|p_{1r}(h)|, \dots, |p_{Mr}(h)|)^\mathrm{T} \quad (10)$$

and

$$M_2(h) = \frac{1}{R' - 1} \sum_{r=1}^{R'} [p_{1r}(h) - p_1(h), \dots, p_{Mr}(h) - p_M(h)] \times [p_{1r}(h) - p_1(h), \dots, p_{Mr}(h) - p_M(h)]^\mathrm{H}, \quad (11)$$

in which $p_m(h) = (1/R') \sum_{r=1}^{R'} p_{mr}(h)$. Here, the sample size $R'$ of simulated measurement should be large enough to guarantee that both moments converge, which can be justified if marginal change of each moment with respect to $R'$ is less than a given threshold. Then, both moments $M_1$ and $M_2$ are interpolated over the source region $\Gamma_1$, such that the values of the first- and second-order moments at every point in $\Gamma_1$ are assigned, denoted by $M_1(\mathbf{r})$ and $M_2(\mathbf{r})$, $\mathbf{r} \in \Gamma_1$.

Finally, the sound source location $\mathbf{r}_s$ is estimated by minimizing the statistical moment error. In the cases of first-order raw and second-order central moments, they are

$$\hat{\mathbf{r}}_s = \arg \min_{\mathbf{r}} \| M_1 - M_1(\mathbf{r}) \|, \quad \mathbf{r} \in \Gamma_1 \quad (12)$$

and

$$\hat{\mathbf{r}}_s = \arg \min_{\mathbf{r}} \| M_2 - M_2(\mathbf{r}) \|, \quad \mathbf{r} \in \Gamma_1, \quad (13)$$

respectively, in which $\| \cdot \|$ represents the Frobenius norm of a vector or a matrix.

## IV. NUMERICAL EXAMPLE

In this section, a 2D numerical example is investigated to illustrate the model and the proposed source estimation method. The sound source is located at $\mathbf{r}_s = (10, 10)$ (Ref. 44) and 21 microphones are placed in the line $x=5$ with same spacing distances [the $m$th microphone's coordinate is $\mathbf{r}'_m = (5, m-1), m=1, \dots, 21$]. Here, the region $\Gamma = \{(x, y): x, y \in [0,20]\}$ is considered, in which the bulk modulus $\kappa(\mathbf{r})$ is assumed to follow a gamma random distribution with mean $\mu_\kappa = 4$ and standard deviation $\sigma_\kappa = 0.6$, and obtained from a Gaussian germ:

$$\kappa(\mathbf{r}) = F_G^{-1}(\xi_{\mathbf{r}}; \mu_{\kappa}, \sigma_{\kappa}^2), \quad (14)$$

where $\xi_{\mathbf{r}}$ is a standard Gaussian random field with covariance function $C(\mathbf{r}_1, \mathbf{r}_2) = e^{(\|\mathbf{r}_1 - \mathbf{r}_2\|)/\ell_c}$ (the correlation length $\ell_c = 10$ m) and $F_G(\cdot; \mu_\kappa, \sigma_\kappa^2)$ is the cumulative function of Gamma distribution with mean $\mu_\kappa$ and variance $\sigma_\kappa^2$. It should be pointed out that Eq. (14) just slightly changes the target correlation function of the random field $C(\mathbf{r}_1, \mathbf{r}_2)$ (cf. Ref. 38). Then, the random sound speed field in the considered region is obtained by $v(\mathbf{r}) = \sqrt{\kappa(\mathbf{r})/\rho}$, in which the density $\rho = 1$ is assumed as a constant. Figure 2(a) shows the

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) (a) A realization of gamma random field $\kappa(\mathbf{r})$ of bulk modulus; (b) wavefront in a random medium at snapshot $t = 6$ s. The point and 21 crosses stand for the sound source and microphones, respectively.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) (a) Five Lagrange polynomials of order 4; (b) GLL points of degree 4.</figcaption></figure>

experimental setup and a realization of the gamma random field $κ(\mathbf{r})$ of bulk modulus. Figure 2(b) displays a wavefront in an (randomly generated) inhomogeneous medium at snapshot $t=6$ s: In this case, the sound speed and wavelength are irregular in different directions due to the inhomogeneity of sound field.

In this work, the sound propagation process in the inhomogeneous medium is carried out via the SEM, in which Lagrange interpolation functions of degree 4–10 are typically used. In this study, Lagrange polynomials of order $n=7$ are employed to perform the numerical simulations of wave propagation. The control points $x_i$ in the $x$-direction ($i \in \{0, 1, 2, \dots, n\}$) are chosen to be the $n$ Gauss Lobatto Legendre (GLL) points, which are the solutions of the following equation:

$$ (1-x^2) \frac{dP_n(x)}{dx} = 0, \qquad (15) $$

in which $[dP_n(x)]/dx$ is the derivative of the Legendre polynomial of order $n$

$$ P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} [(x^2 - 1)^n]. \qquad (16) $$

Using $n$ control points in each direction results in $n^2$ degrees of freedom over each element which are not uniformly spaced. Figure 3(a) depicts the Lagrange polynomials of degree 4, all of which pass through the GLL points, which are the five roots of the equation $(1-x^2)P'_4(x) = 0$, i.e., $x_{GLL} = \{-1, -0.65, 0, +0.65, +1\}$. These points create a mesh, which is shown in Fig. 3(b). It is remarkable that at least five GLL points per minimum wavelength are needed to guarantee the accurate simulation of the SEM.$^{39–41}$ In the present experiment, the average wavelength is $E[v(\mathbf{r})]/f = 2$ m, thus in each direction the considered region must be discretized by at least 50 points. In this work, each unit interval is discretized by eight GLL points (the order of Lagrange polynomial is $n=7$). Therefore, the considered region $\Gamma$ (20 m × 20 m) is discretized by a grid $\Gamma^\#$ having $I = 141 \times 141$ points, which ensures the precision of sound propagation simulation and source localization. Here, a standard Gaussian random field $\xi_\mathbf{r}$ is generated, whose mean $\mu=0$ and variance $\sigma^2=1$ are assumed to be constant. The covariance function $C(\mathbf{r}_1, \mathbf{r}_2) = e^{(\|\mathbf{r}_1 - \mathbf{r}_2\|)/\ell_c}$, where $\mathbf{r}_i \in \Gamma^\#$ for $i=1, 2$ and the correlation length $\ell_c = 10$ m. Figure 4 shows the first 30 eigenvalues of the $141^2 \times 141^2$ covariance matrix $[C(\mathbf{r}_1, \mathbf{r}_2)]_{\mathbf{r}_1, \mathbf{r}_2=1}^I$. This figure justifies that the eigenvalue tends

to zero (from around 8000) within a small number of terms, which implies that most terms in the K L expansion Eq. (7) can be ignored. Here, only the first 50 terms are kept (i.e., $N=50$) such that the normalized global mean square error on random field reconstruction [cf. Eq. (6)] is less than 0.05. Therefore, the random field is described by 50 independent and identically distributed Gaussian random variables.

In this example a Ricker wavelet$^{42}$ with fundamental frequency $f=1$ Hz is used to generate the source signal. Then, $R$ groups of simulated sound pressures are obtained using the SEM and recorded by the $M$ microphones. The recording time is 7.5 s and the sampling frequency is 200 Hz. After the Fourier transform, the measurements at frequency $f=1$ Hz are obtained and denoted by $\mathbf{p} = \{p_{mr}: m=1, \dots, M, r=1, \dots, R\}$. Let $\mathbf{p}(\xi) = [p_1(\xi), \dots, p_M(\xi)]^\mathrm{T}$ denote the random vector of the $M$ measurements. Then, the PDFs, as well as the statistical moments, of $\mathbf{p}(\xi)$ can be computed via the measurement $\mathbf{p}$. Figure 5 shows the PDFs of the real and imaginary parts of the sound pressures from 11 microphones whose $y$-coordinates are 0, 2, 4, ..., 20. Furthermore, the convergence properties of the sample expectation [Eq. (8)] and covariance matrix [Eq. (9)] with respect to the sample size $R$ are shown in Fig. 6. In the left plot, the sample expectations with sample size $R=1, \dots, 500$ from five microphones (y-coordinates are 0, 5, 10, 15, 20) are displayed. The right plot shows the corresponding results for the sample variance. Both plots illustrate that all the sample statistical moments become stable within a

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) The 30 largest eigenvalues of covariance matrix of the random field $\xi_\mathbf{r}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) PDFs of real (left) and imaginary (right) parts of sound pressures from 11 microphones (the y coordinates $y_m$ are 0, 2, 4,..., 20).</figcaption></figure>

sample size of about $R=300$. In this figure, the moments from microphones at $y_m=0$ and $y_m=20$ (correspondingly $y_m=5$ and $y_m=15$) tend to the same values since the microphones are symmetric about the sound source.

Assume that the distance between the sound source and the microphone plane is known, the y-coordinate of the sound source is then estimated. First, the case of sample size $R = R' = 300$ is considered,⁴⁵ the source localization process includes the following three steps:

(1) *Searching the source from a grid with sparse points.* First, the source location is roughly searched from a large region $\Gamma_1 = \{(x, y) : x = 5, y \in [0, 20]\}$. It is averagely discretized by a grid with 21 points, denoted as

$\Gamma_1^\# = \{\mathbf{r}_h = (x_h, y_h) : x_h = 5, y_h = 1, 2, \dots, 21\}$. For each assumed source $\mathbf{r}_h \in \Gamma_1^\#$, the acoustical measurements from the microphones are simulated $R'$ times and the corresponding statistical moments Eqs. (10) and (11) are computed. Then, the source location can be estimated using the minimum error of first- and second-moments, respectively. Figures 7(a) and 7(b) show the errors of both moments between the real measurement and the simulated measurement whose assumed sources are from $\Gamma_1^\#$. Both methods show that the source coordinate is in the interval $y \in [8,12]$.

(2) *Searching the source from a grid with dense points.* Then, the possible region of source can be narrowed: some denser points (but in a smaller region) in the source

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Sample expectations (left) and variances (right) of the measurements from five microphones (the y coordinates $y_m$ are 0, 5, 10, 15, 20), $R$ is the sample size.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) Errors of expectation [subfigures (a), (c), and (e)] and covariance matrix [subfigures (b), (d), and (f)] between the real measurement and the simulated measurement from an assumed source. The y coordinates of the assumed sources are from 21 points in [0,20] [subfigures (a) and (b)], 25 points in [8,12] [subfigures (c) and (d)], and 13 points in [8,12] [sub figures (e) and (f)], respectively. The points represent the errors of the discrete points from the grid and the solid lines stand for the errors of continuous source coordinates obtained by the spline interpolation method. Sample size *R* *R*′ 300. The sound pressures are measured from 21 micro phones with y coordinates 1, 2,..., 20.</figcaption></figure>

plane are chosen: a grid with 25 GLL points whose y-coordinates belong to $y \in [8,12]$ is considered as $\Gamma_1^#$. The corresponding results are shown in Figs. 7(c) and 7(d): the points represent the measurement errors for the assumed sources from $\Gamma_1^#$. In this case, the source location estimates using the first- and second-order moment methods are 10.08 and 9.915 (the exact value is 10), respectively.

(3) *Searching the source from a continuous region.* Finally, both statistical moments are interpolated using the cubic spline interpolation method<sup>43</sup> (which returns a smooth function and is more accurate than linear interpolation), such that the moment errors take values in the continuous interval $y \in [8,12]$ [solid lines in Fig. 7(c) and 7(d)]. The source location is then estimated via Eqs. (12) and (13). In this case, the estimated y-coordinates of the source are, respectively, 10.06 and 10.05, which imply that the estimates are improved by using the interpolation approach.

In order to see how the point density of grid $\Gamma_1^#$ affects the estimation accuracy, $\Gamma_1^#$ is chosen as 13 points whose

y-coordinates belong to $y \in [8,12]$. Above Steps 2 and 3 are then repeated: the corresponding results are shown in Figs. 7(e) and 7(f). In this case, the estimates of source are 10.1 for the expectation error and 10.15 for the covariance matrix error [the minimum points of the solid lines in Figs. 7(e) and 7(f)]; both estimates are less accurate than previously but still acceptable.

Then, a smaller sample size $R = R' = 100$ is considered: the above procedure is repeated and the corresponding results are shown in Fig. 8. In this case, the estimated y-coordinates of source using the first-order statistical moment are 9.915 [from discrete points of $\Gamma_1^#$, the points in subfigure (c)] and 9.94 [from continuous points in $\Gamma_1$, the solid line in subfigure (c)]. Similarly, the estimates using the second-order moment are shown in Fig. 8(d): the estimate is also improved by using the interpolation method, from 10.08 to 10.06. As previously, reducing the number of assumed sources of $\Gamma_1^#$ also decreases the accuracy of estimation. Figures 8(e) and 8(f) show that in the case of $\Gamma_1^#$ being 13 points in $y \in [8,12]$, the source estimates using both moment methods are 11.59 and 10.14, respectively. In particular, under this condition, the expectation error (using both discrete grid and
<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Errors of expectation [subfigures (a), (c), and (e)] and covariance matrix [subfigures (b), (d), and (f)] between the real measurement and the simulated measurement from an assumed source. The y coordinates of the assumed sources are from 21 points in [0,20] [subfigures (a) and (b)], 25 points in [8,12] [subfigures (c) and (d)], and 13 points in [8,12] [sub figures (e) and (f)], respectively. The points represent the errors of the discrete points from the grid and the solid lines stand for the errors of continuous source coordinates obtained by the spline interpolation method. Sample size *R* *R*′ 100. The sound pressures are measured from 21 microphones with y coordinates 1, 2,..., 20.</figcaption></figure>

continuous interval) cannot satisfactorily localize the source. Therefore, either the sample size *R*′ or the density of grid Γ\# must be increased.

Actually, in the case of smaller sample size, the source can still be accurately localized [except the case of Fig. 8(e)]. However, the estimation is less robust; the curve which represents the error is more irregular. This fact can be alternatively justified by the convergence of the statistical moments showed in Fig. 6. Therefore, a reasonable choice of the sample size (a trade-off between estimation robustness and computational cost) may be the one where the concerned statistics become stable, say 300. On the other hand, both Figs. 7 and 8 show that the second-order moment method is more robust than the first-order moment: the error reaches minimum more clearly and has less local fluctuations. This phenomenon may be interpreted as the second moment includes more information (e.g., the relevance between the microphones) than the other.

According to the spatial Nyquist criterion, in order to avoid the spatial aliasing, the spacing distances of the microphones should be smaller than half-wavelength. In the above experiments, the spacing is equal to the average

half-wavelength (1 m), thus the microphone array is adapted to the experiments. In order to investigate the limit of the proposed method, another experiment with fewer microphones is proposed. In this case, only six microphones located at (10, y_m), y_m = 0, 4, 8, 12, 16, 20, are assumed, therefore the spacing is equal to two-wavelength of the sound wave. Let the sample size be *R* = *R*′ = 300, the corresponding computation results are shown in Fig. 9. In this case, the estimated y-coordinates of source with 21 discrete points in Γ\# are 9.96 and 10.05 [minima of the solid lines in subfigures (c) and (d)], respectively. Similarly, the estimates using 13 points in Γ\# are 9.91 and 10.14 for both moments [minima of the solid lines in subfigures (e) and (f)]. It is clear that in this case the estimation results are still precise. However, the proposed method seems less robust (than the results in Fig. 7) since more fluctuations appear in Fig. 9.

Finally, Fig. 10 displays the relative errors of source estimates using both moment methods [lines (a), (c), and (e) for the expectation approach and lines (b), (d), and (f) for the covariance matrix approach] with respect to the number of discrete points in the grid Γ\# (*H* = 7, 13, 19, 25). The results
<figure><img src="image_10.png" /><figcaption>FIG. 9. (Color online) Errors of expectation [subfigures (a), (c), and (e)] and covariance matrix [subfigures (b), (d), and (f)] between the real measurement and the simulated measurement from an assumed source. The y coordinates of the assumed sources are from 21 points in [0,20] [subfigures (a) and (b)], 25 points in [8,12] [subfigures (c) and (d)], and 13 points in [8,12] [sub figures (e) and (f)], respectively. The points represent the errors of the discrete points from the grid and the solid lines stand for the errors of continuous source coordinates obtained by the spline interpolation method. Sample size *R* *R*′ 300. The sound pressures are measured from six microphones with y coordinates 0, 4, 8, 12, 16, 20.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. (Color online) Relative errors of sound source estimates with respect to the point density of grid $\Gamma_1$ using the expectation [(a), (c), and (e)] and covariance matrix [(b), (d), and (f)] methods, for the cases *R* *R*′ 300 and *M* 21 [(a) and (b)], *R* *R*′ 100 and *M* 21 [(c) and (d)], and *R* *R*′ 300 and *M* 6 [(e) and (f)].</figcaption></figure>

for all the above three cases $R = R' = 300$ and $M = 21$ [(a) and (b)], $R = R' = 100$ and $M = 21$ [(c) and (d)], and $R = R' = 300$ and $M = 6$ [(e) and (f)] are presented. This figure shows that the estimation errors of both methods tend to increase as $H$ decreases. All the estimates are acceptable except the one using the expectation error method and low point density of $\Gamma_1^#$. Besides, the covariance matrix error method is generally more accurate than the expectation approach.

V. CONCLUSIONS AND PERSPECTIVES

This paper addresses the problem of sound source localization from the measurements obtained by an array of microphones. The sound propagation medium is assumed to be randomly inhomogeneous, which means that the acoustical field cannot be measured but only its statistical behavior (expectation and covariance function) is known. Instead of an infinite dimensional random field, the medium is simply represented by a small number of uncorrelated random variables via the K L expansion. By sampling from the random field and simulating the sound propagation using the SEM,
the statistical properties (PDF and statistical moments) of
sound pressure at any point in the considered region can be
obtained. Finally, the sound source is localized using a
matched statistical moment method: the source location esti-
mate is decided by minimizing the error of statistical
moment between the real measurement and the simulated
measurement from an assumed source.

A numerical example is presented to illustrate the model
and justify the proposed source localization method. This
experiment indicates that the random sound field can be
effectively replicated by a small number of uncorrelated and
identically distributed random variables using the K L
expansion. The convergence property of statistical moments
is also investigated. Finally, the sound source location is esti-
mated using the expectation and covariance matrix errors of
measurements, respectively. The source localization meth-
ods are proved to be efficient and the sample size needed in
both methods is acceptable.

In general, the wave equation has no analytical solution
in the case of inhomogeneous medium, therefore this work
numerically computes the sound pressure at each point of
space for each sound source, which actually constructs a nu-
merical Green’s function of sound propagation equation. As
a matter of fact, in some special cases of propagation me-
dium (for example, linearly inhomogeneous), a “uniform”
Green’s function (either analytical or numerical) with
respect to all the possible source positions may be con-
structed. Then, the computational cost can be further
reduced. Moreover, real applications could be considered via
the proposed model and source localization method, for
example in underwater and aeroacoustic media, which could
be reasonably assumed as randomly inhomogeneous. These
issues may be conducted in the future works.

ACKNOWLEDGMENTS

This work has been carried out thanks to the support of
the A*MIDEX grant (reference number: ANR-11-IDEX-
0001-02) funded by the French Government “Investissements
d’Avenir” program.

<footnote n="1">¹H. Krim and M. Viberg, “Two decades of array signal processing research,” *IEEE Signal Process. Mag.* **13**(4), 67–94 (1996).</footnote>

<footnote n="2">²B. D. Van Veen and K. M. Buckley, “Beamforming: A versatile approach to spatial filtering,” *IEEE Acoust., Speech Signal Process. Mag.* **5**, 4–24 (1988); available at http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=665.</footnote>

<footnote n="3">³J. Hald, “Basic theory and properties of statistically optimized near field acoustical holography,” *J. Acoust. Soc. Am.* **125**(4), 2105–2120 (2009).</footnote>

<footnote n="4">⁴R. Steiner and J. Hald, “Near field acoustical holography without the errors and limitations caused by the use of spatial DFT,” Int. J. Acoust. Vib. **6**(2), 83–89 (2001); available at http://www.scopus.com/record/display.uri?eid=2 s2.0 0042517094&origin=inward&txGid=0.</footnote>

<footnote n="5">⁵E. G. Williams, *Fourier Acoustic: Sound Radiation and Nearfield Acoustical Holography* (Academic, New York, 1999).</footnote>

<footnote n="6">⁶P. Blomgren, G. Papanicolaou, and H. Zhao, “Super resolution in time reversal acoustics,” *J. Acoust. Soc. Am.* **111**(1), 230–248 (2002).</footnote>

<footnote n="7">⁷P. Druault, R. Marchiano, and P. Sagaut, “Localization of aeroacoustic sound sources in viscous flows by a time reversal method,” *J. Sound Vib.* **332**, 3655–3669 (2013).</footnote>

<footnote n="8">⁸M. Fink, “Time reversed acoustics,” *Sci. Am.* November, 91 97 (1999); avail
able at http://www.scientificamerican.com/article/time reversed acoustics/.</footnote>

<footnote n="9">⁹M. Fink and C. Prada, “Acoustic time reversal mirrors,” *Inverse Probl.* **17**(1), 1–38 (2001).</footnote>

<footnote n="10">¹⁰G. Papanicolaou, K. Solna, and L. Ryzhik, “Statistical stability in time reversal,” *SIAM J. Appl. Math.* **64**(4), 1133–1155 (2004).</footnote>

<footnote n="11">¹¹E. Vergnault, O. Malaspinas, and P. Sagaut, “Noise source identification with the lattice Boltzmann method,” *J. Acoust. Soc. Am.* **133**(3), 1293–1305 (2013).</footnote>

<footnote n="12">¹²S. Finette, “A stochastic representation of environmental uncertainty and its coupling to acoustic wave propagation,” *J. Acoust. Soc. Am.* **120**(5), 2567–2579 (2006).</footnote>

<footnote n="13">¹³J. Ko, D. Lucor, and P. Sagaut, “Effects of base flow uncertainty on Couette flow stability,” *Comput. Fluids* **43**, 82–89 (2011).</footnote>

<footnote n="14">¹⁴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. Am.* **86**(5), 1950–1959 (1989).</footnote>

<footnote n="15">¹⁵E. Shorey, G. De Roeck, and G. Lombaert, “Dealing with uncertainty in model updating for damage assessment: A review,” *Mech. Syst. Signal Process.* **56–57**, 123–149 (2015).</footnote>

<footnote n="16">¹⁶J. Tabrikian and J. L. Krolik, “Robust maximum likelihood source local
ization in an uncertain shallow water waveguide,” *J. Acoust. Soc. Am.*
**101**(1), 241–249 (1997).</footnote>

<footnote n="17">¹⁷J. Tabrikian and J. L. Krolik, “Barankin bounds for source localization in an uncertain ocean environment,” *IEEE Trans. Signal Process.* **47**(11), 2917–2927 (1999).</footnote>

<footnote n="18">¹⁸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.* **66–67**, 159–177 (2016).</footnote>

<footnote n="19">¹⁹M. Shinozuka and G. Deodatis, “Simulation of stochastic processes by spectral representation,” *Appl. Mech. Rev.* **44**(4), 191–204 (1991).</footnote>

<footnote n="20">²⁰R. Cottereau, D. Clouteau, J. P. Vilotte, and R. Madariaga, “Validation of software for 3D propagation of waves in heterogeneous and random media,” *9th US National Congress on Computational Mechanics (USNCCM)*, San Francisco, CA (2007).</footnote>

<footnote n="21">²¹S. Huang and X. Kou, “An extended stochastic response surface method for random field problems,” *Acta Mech. Sin.* **23**, 445–450 (2007).</footnote>

<footnote n="22">²²D. Lucor, J. Meyers, and P. Sagaut, “Sensitivity analysis of large eddy simulations to subgrid scale model parametric uncertainty using polyno- mial chaos,” *J. Fluid Mech.* **585**, 255–279 (2007).</footnote>

<footnote n="23">²³D. Xiu and G. E. Karniadakis, “Modeling uncertainty in flow simulations via generalized polynomial chaos,” *J. Comput. Phys.* **187**(1), 137–167 (2003).</footnote>

<footnote n="24">²⁴W. Betz, I. Papaioannou, and D. Straub, “Numerical methods for the dis cretization of random fields by means of the Karhunen Loève expansion,” *Comput. Methods Appl. Mech. Eng.* **271**, 109–129 (2014).</footnote>

<footnote n="25">²⁵K. K. Phoon, S. P. Huang, and S. T. Quek, “Simulation of second order processes using Karhunen Loeve expansion,” *Comput. Struct.* **80**, 1049–1060 (2002).</footnote>

<footnote n="26">²⁶M. V. De Hoop and K. Solna, “Estimating a Green’s function from “field field” correlations in a random medium,” *SIAM J. Appl. Math.* **69**(4), 909–932 (2009).</footnote>

<footnote n="27">²⁷Y. L. Li, C. H. Liu, and S. J. Franke, “Three dimensional Green’s function for wave propagation in a linearly inhomogeneous medium the exact analytic solution,” *J. Acoust. Soc. Am.* **87**(6), 2285–2291 (1996).</footnote>

<footnote n="28">²⁸V. Vavrycuk, “Acoustic and elastodynamic 3D Green’s functions for isotropic media with a weak velocity gradient,” *Wave Motion* **31**, 223–236 (2000).</footnote>

<footnote n="29">²⁹J. F. Semblat, “Modeling seismic wave propagation and amplification in 1D/2D/3D linear and nonlinear unbounded media,” *Int. J. Geomech.* **11**(6), 440–448 (2010).</footnote>

<footnote n="30">³⁰A. T. Patera, “A spectral element method for fluid dynamics: Laminar flow in a channel expansion,” *J. Comput. Phys.* **54**(3), 468–488 (1984).</footnote>

<footnote n="31">³¹Y. Maday and A. T. Patera, “Spectral element methods for the Navier Stokes equations,” in *State of the Art Survey in Computational Fluid Mechanics*, edited by A. K. Noor and J. T. Oden (ASME, New York, 1988), Chap. 3.</footnote>

<footnote n="32">³²E. Faccioli, F. Maggio, R. Paolucci, and A. Quarteroni, “2D and 3D elastic wave propagation by a pseudo spectral domain decomposition method,” *J. Seismol.* **1**(3), 237–251 (1988).</footnote>

<footnote n="33">³³D. Komatitsch and J. P. Vilotte, “The spectral element method: An effi-
cient tool to simulate the seismic response of 2D and 3D geological
structures,” Bull. Seismol. Soc. Am. **88**(2), 368–392 (1998); available at
http://bssa.geoscienceworld.org/content/88/2/368.short</footnote>

<footnote n="34">³⁴R. Paolucci, E. Faccioli, and F. Maggio, “3D response analysis of an instrumented hill at Matsuzaki, Japan, by a spectral method,” *J. Seismol.* **3**(2), 191–209 (1999).</footnote>

<footnote n="35">³⁵E. Chaljub, D. Komatitsch, J. P. Vilotte, Y. Capdeville, B. Valette, and G. Festa, “Spectral element analysis in seismology,” *Adv. Geophys.* **48**, 365–419 (2007).</footnote>
<footnote n="36"><sup>36</sup>N. E. Horlin, “3D hierarchical hp FEM applied to elasto acoustic model ling of layered porous media,” *J. Sound Vib.* **285**(1), 341–363 (2005).</footnote>

<footnote n="37"><sup>37</sup>O. P. Maitre and O. M. Knio, *Spectral Method for Uncertainty Quantification With Applications to Computational Fluid Dynamic* (Springer, New York, 2010).</footnote>

<footnote n="38"><sup>38</sup>M. Grigoriu, “Simulation of stationary non Gaussian translation pro-esses,” *J. Eng. Mech.* **124**(2), 121–126 (1998).</footnote>

<footnote n="39"><sup>39</sup>P. Cristini and D. Komatitsch, “Some illustrative examples of the use of a spectral element method in ocean acoustics,” *J. Acoust. Soc. Am.* **131**(3), EL229 EL235 (2012).</footnote>

<footnote n="40"><sup>40</sup>G. Seriani and E. Priolo, “Spectral element method for acoustic wave simula-tion in heterogeneous media,” *Finite Elem. Anal. Des.* **16**(3), 337–348 (1994).</footnote>

<footnote n="41"><sup>41</sup>J. Tromp, D. Komattisch, and Q. Liu, “Spectral element and adjoint methods in seismology,” Commun. Comput. Phys. **3**(1), 1–32 (2008); available at http:// resolver.caltech.edu/ CaltechAUTHORS: TROccp08.</footnote>

<footnote n="42"><sup>42</sup>J. F. Semblat, L. Lenti, and A. Gandomzadeh, “A simple multi directional absorbing layer method to simulate elastic wave propagation in unbounded domains,” *Int. J. Numer. Methods Eng.* **85**, 1543–1563 (2011).</footnote>

<footnote n="43"><sup>43</sup>O. Dubrule, “Comparing splines and kriging,” *Comput. Geosci.* **10**, 327–338 (1984).</footnote>

<footnote n="44"><sup>44</sup>In this section, the spatial coordinates are in meter (m), except where noted.</footnote>

<footnote n="45"><sup>45</sup>In the experiments of this section, the sample sizes *R* and *R*′ of real and simulated measurements are assumed to be the same, but this is not an essential condition.</footnote>