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 . 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 . 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 
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 
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 distributionthat 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 . 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 constantis
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  that this well-posedness requires the continuity of the level set map. In , this point has been discussed and is given in the following proposition
 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 . 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 
are the eigenvalues and
are the eigenfunctions of the compact and self-adjoint integral operatordefined 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 .
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
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
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.
For fixed , acts as a bounded linear functional on and therefore is a real valued Gaussian random variable, which implies . 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. ∎
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.
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
5 Numerical test
We characterize the posterior distribution by means of sampling with the preconditioned Crank-Nicolson (pCN) MCMC method developed in . In concrete, this algorithm is implemented in the following iteration process:
|For , initialize|
|For , propose||,||,|
|where .||where .|
|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 . 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.
. 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.
-  A. Aghasi, M. Kilmer, and E. Miller, Parametric level set methods for inverse problems, SIAM J. Imaging Sciences, 4(2), (2011), 618-650.
-  A. Alzaalig, G. Hu, X. Liu and J. Sun, Fast acoustic source imaging using multi-frequency sparse data, arXiv:1712.02654.
-  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.
-  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.
-  M. Burger, A level set method for inverse problems, Inverse problems, 17(5): 1327-1355, 2001.
-  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.
-  N. Chada, M. Iglesias, L. Roininen and A. Stuart, Parameterizations for ensemble Kalman inversion, Inverse Problems 34 (2018) 055009 (31pp).
-  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.
-  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.
-  O. Dorn and D. Lesselier, Level set methods for inverse scattering, Inverse Problems 22 (2006) R67–R131.
-  M. Dunlop, M. Iglesias and A. Stuart, Hierarchical Bayesian level set inversion, Statistics and Computing, November 2017, Volume 27, Issue 6, 1555-1584.
-  M. Dunlop and A. Stuart, The Bayesian formulation of EIT: analysis and algorithms, Inverse Problems and Imaging, 10(4), 2016, 1007–1036.
-  M. Eller and N. Valdivia, Acoustic source identification using multiple frequency information, Inverse Problems 25, 115005(20pp), 2009.
-  R. Griesmaier, Multi-frequency orthogonality sampling for inverse obstacle scattering problems, Inverse Problems, 27 (2011), 085005.
-  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
-  M. Iglesias, Y. Lu and A. Stuart, A Bayesian Level Set Method for Geometric Inverse Problems, Interfaces and Free Boundaries 18 (2016), 181-217.
-  M. Iglesias, K. Law and A. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems 29 (2013) 045001 (20pp).
-  M. Iglesias, A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems, Inverse Problems 32 (2016) 025002 (45pp).
-  M. Iglesias, K. Lin, A. Stuart, Well-posed Bayesian geometric inverse problems arising in subsurface flow, Inverse Problems, 30 2014, 114001.
-  K. Ito, K. Kunisch and Z. Li, Level-set function approach to an inverse interface problem, Inverse Problems 17 (2001) 1225–1242.
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.
-  J. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer Verlag, London, 2005.
-  J. Kaipio, V. Kolehmainen, M. Vauhkonen and E. Somersalo, Inverse problems with structural prior information, Inverse Problems, 15(3), 1999, 713–729.
-  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.
-  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.
-  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.
-  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.
-  X. Liu and J. Sun, Reconstruction of Neumann eigenvalues and support of sound hard obstacles, Inverse Problems, 30 (2014), 065011.
R. Lorentzen, K. Flornes and G. Nævdal, History matching channelized reservoirs using the ensemble Kalman filter. SPE Journal 17 (2012), 137–151.
-  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.
-  S. Osher and R. Fedkiw, The Level Set Method and Dynamic Implicit Surfaces. Springer, New York (2003). Zbl pre01873119
-  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
-  J. Ping and D. Zhang, History matching of channelized reservoirs with vector-based level-set. SPE Journal 19 (2014), 514–529.
-  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
-  R. Potthast, A study on orthogonality sampling, Inverse Problems, 26 (2010), 074015.
C. Rasmussen and C. Williams, ”Gaussian Processes for Machine Learning”, Adaptive Computation and Machine Learning, The MIT Press, 2006.
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.
-  G. Santin and R. Schaback, Approximation of eigenfunctions in kernel-based spaces, Adv Comput Math 42, (2016), 973–993.
-  J. Sun, An eigenvalue method using multiple frequency data for inverse scattering problems, Inverse Problems, 28 (2012), 025012.
-  A. Stuart, Inverse problems: A Bayesian perspective, Acta Numerica (2010), 451–559.
-  F. Santosa, A level-set approach for inverse problems involving obstacles, ESAIM: Control, Optimization and Calculus of Variations, 1: 17-33, 1996.
-  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.
-  A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, 2005.
-  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.
-  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.
-  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.
-  R. Zhang and J. Sun, Efficient finite element method for grating profile reconstruction, J. Comput. Phys., 302 (2015), 405-419.