A parametric Bayesian level set approach for acoustic source identification using multiple frequency information

07/19/2019 ∙ by Zhiliang Deng, et al. ∙ Southwest Jiaotong University NetEase, Inc 0

The reconstruction of the unknown acoustic source is studied using the noisy multiple frequency data on a remote closed surface. Assume that the unknown source is coded in a spatial dependent piecewise constant function, whose support set is the target to be determined. In this setting, the unknown source can be formalized by a level set function. The function is explored with Bayesian level set approach. To reduce the infinite dimensional problem to finite dimension, we parameterize the level set function by the radial basis expansion. The well-posedness of the posterior distribution is proven. The posterior samples are generated according to the Metropolis-Hastings algorithm and the sample mean is used to approximate the unknown. Several shapes are tested to verify the effectiveness of the proposed algorithm. These numerical results show that the proposed algorithm is feasible and competitive with the Matérn random field for the acoustic source problem.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Sound source detection has been the active research topics in the field of acoustic scattering problems due to the extensive applications to various fields such as military surveillance to detect flying objects or armored vehicles and a humanoid robot auditory system. A main object to be determined for this kind of problems is the location and shape of the sound source. These information are coded in a spatial dependent function with compact support. For reconstructing this function, the multiple frequency sound wave data in a remote close surface are usually collected. The observed data is related to the unknown function by a partial differential equation. The task of such source discovery is an inverse problem. It is usually ill-posed, i.e., the solution is highly sensitive to changes in the observed data. A standard way of solving inverse problems for seeking the source function is to minimize some functional related with the collected data and simulated data computed with a candidate source function. Some iterative procedure is typically applied to achieve the minimizer of the functional. To cope with the numerical instability or make the computation more stable, a regularized term is imposed on the minimization function, named the regularization method.

In this paper, we consider using the Bayesian level set technique to reconstruct the unknown acoustic source function with compact support. The level set approach is proposed in [31, 32] to track the wave front interface. Compared with the classical curve parameterization method, the level set approach allows for topological changes to be detected during the course of algorithms. This advantage enables it to become the popular technique to deal with interface problems. In inverse problems, it has also been discussed extensively [5, 10, 20, 41, 42]

. It can be seen that these studies cope with the shape reconstruction problem in the classical framework, rarely analyze the uncertainty of the unknown variables. The uncertainty is often characterized by some statistical techniques. And moreover the characterization provides the estimation for the unknown variables and their reliability information. To quantify the uncertainty, the Bayesian level set algorithm is studied in

[11, 16]. One can refer to [29, 30, 33, 42] for some related applications of this algorithm. In Bayesian level set inversion, the reconstruction target is the posterior distribution, which provides us a basis for quantifying the uncertainty. In this process, the prior distribution needs to be treated firstly, actually before the data is acquired. We use the Whittle-Matérn random field as the prior information of the level set function. Such random field has extensive applications and is studied by many authors [26, 36, 37]. In [7, 11], it is used to characterize the prior of the level set function. As seen in these references, the simulation of the random field can be obtained by solving some related stochastic differential equation or utilizing the Karhunen-Loève expansion by the eigen-system of the covariance function. The stochastic differential equation can be solved by the finite difference or the finite finite approach [26, 36, 37]. In regular domain, the eigen-system can be in general obtained directly [7, 11]. However, in non-regular domain, it is still necessary to use some numerical approximation algorithms to obtain the approximated eigen-system [21]. It can be seen that the prior discretization, whatever by the stochastic differential equation or the Karhunen-Loève expansion, will induce a high dimensional unknown variable. This may cause huge computation burden in Bayesian inversion process. To reduce the computation complexity, some acceleration algorithms have been explored for statistical inversion. It is well-known that the surrogate model algorithms can reduce the numerical complexity effectively [9]. The surrogate model methods try to use the polynomial chaos expansion to replace the forward model or the likelihood function. And in doing so, we just need to evaluate the polynomials instead of the forward model or likelihood function, which significantly reduces the numerical burden in Bayesian inversion.

In our problem, we focus on using the radial basis function (RBF) expansion to parameterize the level set function. The parametric level set approach is considered in the classical frame

[1, 27, 34]. To the authors’ knowledge, it still has not been studied in statistical inversion. We compare this numerical effects for using the Whittle-Matérn random field prior with the RBF expansion prior. And the well-posdeness of the posterior distribution based on the RBF expansion prior is explained.

The remainder is organized in the following: In Section 2, we depict the inverse acoustic wave source reconstruction problem. In Section 3, the Bayesian level set approach based on the Whittle-Matérn random field prior is discussed. In Section 4, we propose the RBF parametric level set algorithm to solve the reconstruction problem. We give some numerical tests to verify the proposed algorithms in Section 5.

2 Acoustic source problem

The pressure of time-harmonic wave radiated from source with density is modeled by [13]


where is the wave number, is the radial frequency, is the speed of sound, is the spatial dimension and (2.2) is the Sommerfeld radiation condition, is a given function and is the unknown with compact support to be determined. The multiple frequency observed data for the acoustic wave field is collected in a remote closed surface . Suppose the acoustic source domain is compactly contained in the region bounded by the closed surface . For simplicity, we consider 2-dimensional case, i.e., and let . According to the potential representation, we have [13]


where is the cylindrical Hankel function. The data are collected on finite points for finite wave numbers . Assume that the noise is additive Gaussian, ,

is the noise mean variance. By (

2.3), the system can be written as


where depends on wave number and is the noise. Denote and . In this paper, we assume that the function is known a priori to have the form


here denotes the indicator function of subset , are subsets of such that and (), the are known positive constants. In this setting, the regions determine the unknown function and therefore become the primary unknowns. Some papers have discussed the acoustic source reconstruction problems [2, 3, 13, 45]. From these existed results, one can see that the domain can be determined uniquely by the multiple frequency data.

3 Level set Bayesian inversion with Whittle-Matérn random field prior

As stated above, instead of reconstructing the acoustic source, we solve the recovery problem of interfaces or regions . It can be devised as an interface tracking problem. A simple and versatile approach for this kind of problems is the level set technique. In this method, the interface is characterized by a so-called level set function. Actually, the interface is represented as the contour curve at some fixed value of the level set function, which is a higher dimensional function. For our problem, the level set function is defined by the following way [10, 11, 16]:


where are constants. By this characterization, the source function is determined. Define the level set map


Since there exist many different determine the same , the level set map is not an injection. However, if the function is fixed, the function is uniquely specified. With (2.4) and (3.2), we have the following operator equation


We deal with the equation in the frame of Bayesian inversion. Bayesian perspective views the involved quantities

as random variables. Denote the random variable

, given

. The solution of Bayesian reconstruction is the posterior probability distribution

that expresses our beliefs regarding how likely the different parameter values are. The posterior density is denoted by . Bayes’ formula shows that


where is the likelihood function, the prior density and the constant. Since the noise is additive Gaussian , the likelihood function is given by


A frequently used prior for the unknown is the Whittle-Matérn random field . The covariance function is given by


where is the smoothness parameter, is the modified Bessel function of the second kind of order and is the length-scale parameter. The integer value of determines the mean square differentiability of the underlying process, which matters for predictions made using such a model. This distribution is discussed by many authors [26, 36, 37]. A method for explicit, and computationally efficient, continuous Markov representations of Gaussian Matérn fields is derived [26]. The method is based on the fact that a Gaussian Matérn field on can be viewed as a solution to the stochastic partial differential equation (SPDE)



is the Gaussian white noise and the constant


Here is the variance of the stationary field. The operator

is a pseudo-differential operator defined by its Fourier transform. When

, we can solve (3.7) by the finite element method with suitable boundary condition. The stochastic weak solution of the SPDE (3.7) is found by requiring that


where . We use the linear element to solve the weak formulation with Dirichlet boundary condition, i.e., .

For Bayesian inversion, the well-posedness of the posterior distribution can be established by some assumptions on the forward operator and the prior distribution on [16, 40]. It can be seen from [16] that this well-posedness requires the continuity of the level set map. In [16], this point has been discussed and is given in the following proposition

Proposition 3.1.

[16] For and , the level set map is continuous at if and only if for all . Here denotes the Lebesgue measure and is the -th level set .

4 Parametric level set based on RBF expansion

Bayesian inversion requires the numerical evaluation of the forward problem repeatedly for a large number of samples. An appropriate parameterization of the unknown function can significantly reduce the dimensionality of an inverse problem. A low order model based on the radial basis expansion for the level set function is discussed in [1, 27, 34], where the level set function is represented by the radial basis function expansion. This representation provides flexibility in terms of its ability to characterize shapes of varying degrees with varying degrees of complexity as measured specifically by the curvature of the boundary [1]. We represent the unknown level set function using the radial basis expansion


where are the random coefficients, the radial basis function and called the RBF centers. Here, we take the Gaussian radial basis function


where are the Gaussian widths at centers . The kernel is Mercer kernel and therefore has the following representation [38]



are the eigenvalues and

are the eigenfunctions of the compact and self-adjoint integral operator

defined by


Moreover, the level set function can be written as


Denote . We assume that the components are independent and identically distributed random variables. Denote the probability space by . Define the map according to (4.1)


The measure space is the push-forward of under . The measure is well-defined since the RBF is continuous. Actually, the new measure is defined by


With the expansion (4.1), we have the recovery problem of the coefficient


This likelihood function is given for (4.8)


Therefore, the posterior density satisfies


where is the prior density for .

Theorem 4.1.

If is -dimensional Gaussian, i.e., , then is the Gaussian measure . The covariance function of this random field is of the form


We can get the covariance function of the random field

By virtue of the parallelogram law, it follows that


where .

Next, we discuss that the RBF expansion of the level set function can meet the condition in Proposition 3.1, i.e., the measure of the level set vanishes almost surely on both and . Define the functional by

Lemma 4.2.

For any , is -measurable.


It suffices to verify for any , the set . Since is non-negative, for , we know that and hence measurable. For , it can be obtained that is open. In fact, if it is not open, then there exists , such that is not an interior point. That means that for every , there exists , even if it holds that , we still have . Therefore we have . Moreover, by the construction of and the triangle inequality, we have . Noting that

and that is decreasing, we can conclude that

This is a contradiction. So is open for . Then we know that is -measurable. ∎

Introduce the random level set


By Lemma 4.2, it follows that is a random variable on both and . Next we demonstrate that vanishes almost surely on these two measure spaces.

Theorem 4.3.

Let the coefficient vector

in (4.1) be drawn from the

-dimensional Gaussian distribution. The level set function is given by (

4.1). The random level set of is defined by (4.13). Then

  1. -almost surely;

  2. -almost surely.


For fixed , acts as a bounded linear functional on and therefore is a real valued Gaussian random variable, which implies [16]. Moreover, it follows that


Since , it follows that , -almost surely. Next, we show that -almost surely. Set

Since is open for any , the set is Borel set and measurable. Therefore, we have

The above equality shows that -almost surely. ∎

Remark 1.

According to Proposition 3.1, the level set map is continuous at the radial basis expansion . With this point, we can obtain the well-posedness of the posterior distribution from the theories in [16].

The well-posedness of the posterior distribution is characterized in the sense of the Hellinger metric. Recall that the Hellinger distance between 􏱄 and 􏱄 is defined as

for any measure 􏱅 with respect to which 􏱄 and are absolutely continuous.

Theorem 4.4.

Denote . The prior for is the -dimensional Gaussian. The posterior distribution is locally Lipschitz with respect to in the Hellinger distance: for all , with , there exists a constant such that



It is obvious that and . In fact, we can conclude that are strictly positive. For instance, for , we have that

By , we know that

Therefore, it follows that


Since is -dimensional Gaussian, the unit ball has strictly positive measure, i.e., . This shows that the posterior distribution is well-defined probability distribution. By the definition of Hellinger distance, we have that


For , by (4.16) we have


For , it follows that


Combining (4.17) and (4.18), we get the estimation (4.15). ∎

5 Numerical test

We characterize the posterior distribution by means of sampling with the preconditioned Crank-Nicolson (pCN) MCMC method developed in [40]. In concrete, this algorithm is implemented in the following iteration process:

For , initialize
For , propose , ,
where . where .
Compute probability
Set , if rand, , if rand,
, otherwise. , otherwise.
Return to the second line.

In the numerical tests, we consider some parameter settings for the model problem given in [13]. Assume that the frequencies vary from Hz to kHz and ms. We suppose that the unknown region is contained in a disk domain and the data receiver is located on . The problem geometry is displayed in Fig. 1. The disk domain is partitioned into triangle element. The forward operator is evaluated by using the Gaussian integral in each triangle element in the domain . The Whittle-Matérn random field is generated by (3.8) in the same mesh. The radial basis centers are taken as the vertices of the triangle elements. In addition, we also use center points in the RBF expansion. The centers are displayed in Fig. 2. We first show several samples drawn according to the Whittle-Matérn random field and the RBF expansion in Fig. 3 with and . The first and third rows in Fig. 3 display the samples for the level set function and the second and fourth rows give the following function

We display some numerical reconstructions for the case of the source with single medium in Fig. 4. The noise level . The prior parameters are taken as in Fig. 3. In the MH algorithm, parameter is taken as . In these numerical examples, we generate samples from the posterior distribution. Some samples in the burn-in process being excluded, the sample means are used as the approximations of the true domains. The first two rows give the single domain case and the last row displays two separate domains.

In Fig. 5, we display the numerical reconstructions for the sources with two different mediums for . In these examples, samples are drawn.

We use center points in the RBF expansion and show the numerical reconstructions in Fig. 6 with iteration steps and .

These numerical reconstructions show that the effectiveness of the MCMC. From Fig. 5, we see that the proposed approach can cope well with some complex domains, which provides us conveniences in many cases. When the center points are reduced to , Fig. 6 shows the RBF expansion also works well. This means that by the RBF parameterization, we can reduce the dimensional number of the unknown variable effectively.

We plot the error functions versus the accepted samples number in Fig. 7 corresponding to Fig. 4-6

. It can be seen that the errors decay at the beginning and subsequently fluctuate within a narrow range. This reflects that the Markov chains become stable after some iterations.

Fig. 1: The problem geometry.
Fig. 2: The centers of the RBF function, left: 2161 center points; right: 25 center points.
Fig. 3: Lines 1, 2: Whittle Matérn random field prior samples and the corresponding ; Lines 3, 4: RBF expansion samples and the corresponding .
(i) Lshape
(ii) Single disk
(iii) Two disks with the same medium value
Fig. 4: The numerical reconstructions for . Left: the true domain; Middle: the reconstructions using Matérn field prior; Right: the reconstructions using the RBF expansion.
(i) Two disks with different medium values
(ii) Two different domains with different medium values
(iii) Embedded domain with different medium values
Fig. 5: The numerical reconstructions with for two kinds of mediums. Left: the true domain; Middle: the reconstructions using Matérn field prior; Right: the reconstructions using the RBF expansion.
(i) Lshape
(ii) Single disk
(iii) Two disks with the same medium value
Fig. 6: The numerical reconstructions with using centers in the RBF expansion for single medium case.
(i) Corresponding to Fig. 4
(ii) Corresponding to Fig. 5
(iii) Corresponding to Fig. 6
Fig. 7: The corresponding model data misfit functions for the accepted samples.


  • [1] A. Aghasi, M. Kilmer, and E. Miller, Parametric level set methods for inverse problems, SIAM J. Imaging Sciences, 4(2), (2011), 618-650.
  • [2] A. Alzaalig, G. Hu, X. Liu and J. Sun, Fast acoustic source imaging using multi-frequency sparse data, arXiv:1712.02654.
  • [3] G. Bao, S. Lu, W. Rundell and B. Xu, A recursive algorithm for multi-frequency acoustic inverse source problems, SIAM J. Numer. Anal., 53 (2015), 1608-1628.
  • [4] A. Beskos, A. Jasra, E. Muzaffer, A. Stuart, Sequential Monte Carlo methods for Bayesian elliptic inverse problems, Stat Comput, 2015, DOI 10.1007/s11222-015-9556-7.
  • [5] M. Burger, A level set method for inverse problems, Inverse problems, 17(5): 1327-1355, 2001.
  • [6] T. Bui-Thanh and O. Ghattas, An analysis of infinite dimensional Bayesian inverse shape acoustic scattering and its numerical approximation, SIAM/ASA Journal on Uncertainty Quantification, 2, 2014, 203-222.
  • [7] N. Chada, M. Iglesias, L. Roininen and A. Stuart, Parameterizations for ensemble Kalman inversion, Inverse Problems 34 (2018) 055009 (31pp).
  • [8] S. Cotter, M. Dashti, J. Robinson and A.M. Stuart, Bayesian inverse problems for functions and applications to fluid mechanics. Inverse Problems 25, 2009, 115008.
  • [9] Z. Deng and X. Yang, Convergence of spectral likelihood approximation based on q-Hermite polynomial in Bayesian inversion, Proceedings of the American Mathematical Society, 2019.
  • [10] O. Dorn and D. Lesselier, Level set methods for inverse scattering, Inverse Problems 22 (2006) R67–R131.
  • [11] M. Dunlop, M. Iglesias and A. Stuart, Hierarchical Bayesian level set inversion, Statistics and Computing, November 2017, Volume 27, Issue 6, 1555-1584.
  • [12] M. Dunlop and A. Stuart, The Bayesian formulation of EIT: analysis and algorithms, Inverse Problems and Imaging, 10(4), 2016, 1007–1036.
  • [13] M. Eller and N. Valdivia, Acoustic source identification using multiple frequency information, Inverse Problems 25, 115005(20pp), 2009.
  • [14] R. Griesmaier, Multi-frequency orthogonality sampling for inverse obstacle scattering problems, Inverse Problems, 27 (2011), 085005.
  • [15] R. Griesmaier and C. Schmiedecke. A factorization method for multi-frequency inverse source problems with spars far field measurements, SIAM J. Imag. Sci., (2017), to appear
  • [16] M. Iglesias, Y. Lu and A. Stuart, A Bayesian Level Set Method for Geometric Inverse Problems, Interfaces and Free Boundaries 18 (2016), 181-217.
  • [17] M. Iglesias, K. Law and A. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems 29 (2013) 045001 (20pp).
  • [18] M. Iglesias, A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems, Inverse Problems 32 (2016) 025002 (45pp).
  • [19] M. Iglesias, K. Lin, A. Stuart, Well-posed Bayesian geometric inverse problems arising in subsurface flow, Inverse Problems, 30 2014, 114001.
  • [20] K. Ito, K. Kunisch and Z. Li, Level-set function approach to an inverse interface problem, Inverse Problems 17 (2001) 1225–1242.
  • [21]

    L. Jiang and N. Ou, Bayesian inference using intermediate distribution based on coarse multiscale model for time fractional diffusion equation, Multiscale Modeling and Simulation 16 ( 2018), 327-355.

  • [22] J. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer Verlag, London, 2005.
  • [23] J. Kaipio, V. Kolehmainen, M. Vauhkonen and E. Somersalo, Inverse problems with structural prior information, Inverse Problems, 15(3), 1999, 713–729.
  • [24] J. Kaipio, V. Kolehmainen, E. Somersalo and M. Vauhkonen, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse Problems 16(5), 2000, 1487-1522.
  • [25] J. van der Kruk, C. Wapenaar, J. Fokkema and P. van den Berg, Improved three-dimensional image reconstruction technique for multi-component ground penetrating radar data, Subsurface Sensing Technologies and Applications, 4(1), 61-99, 2003.
  • [26] F. Lindgren and H. Rue and J. Lindstrm, An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approch, J. R. Statist. Soc. B, 73, (2011), 423-498.
  • [27] D. Liu, A. Khambampati and J. Du, A parametric level set method for electrical impedance tomography, IEEE Transactions on medical imaging 37(2), (2018), 451-460.
  • [28] X. Liu and J. Sun, Reconstruction of Neumann eigenvalues and support of sound hard obstacles, Inverse Problems, 30 (2014), 065011.
  • [29]

    R. Lorentzen, K. Flornes and G. Nævdal, History matching channelized reservoirs using the ensemble Kalman filter. SPE Journal 17 (2012), 137–151.

  • [30] R. Lorentzen, G. Nævdal, and A. Shafieirad, Estimating facies fields by use of the ensemble Kalman filter and distance functions – applied to shallow-marine environments. SPE Journal 3 (2012), 146–158.
  • [31] S. Osher and R. Fedkiw, The Level Set Method and Dynamic Implicit Surfaces. Springer, New York (2003). Zbl pre01873119
  • [32] S. Osher and A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. 79 (1988), 12–49. Zbl 0659.65132 MR 89h:80012
  • [33] J. Ping and D. Zhang, History matching of channelized reservoirs with vector-based level-set. SPE Journal 19 (2014), 514–529.
  • [34] G. Pingen, M. Waidmann, A. Evgrafov and K. Maute, A parametric level-set approach for topology optimization of flow domains, Struct Multidisc Optim 41 (2010), 117–131
  • [35] R. Potthast, A study on orthogonality sampling, Inverse Problems, 26 (2010), 074015.
  • [36]

    C. Rasmussen and C. Williams, ”Gaussian Processes for Machine Learning”, Adaptive Computation and Machine Learning, The MIT Press, 2006.

  • [37]

    L. Roininen, J. Huttunen and S. Lasanen, Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography, Inverse Problems and Imaging, 8(2), 2014, 561–586.

  • [38] G. Santin and R. Schaback, Approximation of eigenfunctions in kernel-based spaces, Adv Comput Math 42, (2016), 973–993.
  • [39] J. Sun, An eigenvalue method using multiple frequency data for inverse scattering problems, Inverse Problems, 28 (2012), 025012.
  • [40] A. Stuart, Inverse problems: A Bayesian perspective, Acta Numerica (2010), 451–559.
  • [41] F. Santosa, A level-set approach for inverse problems involving obstacles, ESAIM: Control, Optimization and Calculus of Variations, 1: 17-33, 1996.
  • [42] X. Tai and T. Chan, A survey on multiple level set methods with applications for identifying piecewise constant functions, International Journal of Numerical Analysis and Modeling 1 (2004), 25–47.
  • [43] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, 2005.
  • [44] Z. Wang, J. Bardsley, A. Solonen, T. Cui and Y. Marzouk, Bayesian inverse problems with priors: a randomize-then-optimize approch, SIAM J. Sci. Comput., 39(5), 2017, S140-S166.
  • [45] X. Wang, Y. Guo, D. Zhang and H. Liu, Fourier method for recovering acoustic sources from multi-frequency far-field data, Inverse Problems, 33 (2017), 035001.
  • [46] X. Yang and Z. Deng, An ensemble Kalman filter approach based on operator splitting for solving nonlinear Hammerstein type ill-posed operator equations, Modern Physics Letters B, 32(28), (2018), 1850335.
  • [47] R. Zhang and J. Sun, Efficient finite element method for grating profile reconstruction, J. Comput. Phys., 302 (2015), 405-419.