An ensemble Kalman filter approach based on level set parameterization for acoustic source identification using multiple frequency information

TL;DR

This paper introduces a novel ensemble Kalman filter method combined with level set techniques for reconstructing unknown acoustic sources from noisy multi-frequency data, demonstrating effective numerical results.

Contribution

The paper presents a new statistical inversion algorithm integrating ensemble Kalman filter with level set methods for acoustic source support determination.

Findings

Effective reconstruction of acoustic sources from noisy data

Successful numerical examples demonstrating method accuracy

Supports multi-frequency data integration

Abstract

The spatial dependent unknown acoustic source is reconstructed according noisy multiple frequency data on a remote closed surface. Assume that the unknown function is supported on a bounded domain. To determine the support, we present a statistical inversion algorithm, which combines the ensemble Kalman filter approach with level set technique. Several numerical examples show that the proposed method give good numerical reconstruction.

Equations

$$(△+k^2)u(x,k)=\mathcal{S}(k,x),x\in {\mathbb{R}}^d,$$

$$\displaystyle\lim_{r\rightarrow\infty}r^{\frac{d-1}{2}}{\big{\{}}\frac{\partial u}{\partial r}-iku{\big{\}}}=0,\,\,r=|x|,$$

$$\mathcal{S}(k,x)=A(k)f(x),\text{supp}f\subset D.$$

$$u(k,x)$$

$$=A(k){\int }_{D}f(y)H_0^{(1)}(k|x-y|)dy,$$

$${\mathcal{H}}_{k_m}f+{\eta }_m=[u(x_j,k_m){]}_{j=1}^N:=b_m\text{for}m=1,2,\cdots M,$$

$$f(x)=\sum _{l=1}^n{w}_l{\mathbb{I}}_{D_l}(x);$$

$$D_l=\{x\in D\mid c_{l-1}\le \phi (x)<c_l\},$$

$$\mathcal{G}(\phi )(x)=f(x)=\sum _{l=1}^n{w}_l{\mathbb{I}}_{D_l}(x).$$

$${\mathcal{H}}_{k_m}\circ \mathcal{G}(\phi )=[u(k_m,x_j){]}_{j=1}^N\text{for}m=1,2,\cdots M.$$

$$\mathcal{K}(\phi )+\eta =b,$$

$$\partial D(t)=\{x:\phi (x,t)=0\}.$$

$$\frac{\partial \phi }{\partial t}=-v\cdot \nabla \phi =-v|\nabla \phi |\cdot \frac{\nabla \phi }{|\nabla \phi |}:=-\mathfrak{v}|\nabla \phi |,$$

$$\mathfrak{v}=J(f{)}^T(b-\mathcal{H}f),$$

$$\frac{\partial \phi }{\partial t}=J(f{)}^T(\mathcal{H}f-b)|\nabla \phi |,$$

$$\phi (x,0)={\phi }_0.$$

$${\phi }_{n+1}=\mathcal{M}({\phi }_n)$$

$${\phi }_{n+1}=\mathcal{M}({\phi }_n),n=0,1,\cdots ,N_t$$

$${\phi }_0\sim {\mu }_0,$$

$$\mathcal{K}({\phi }_{n+1})+{\eta }^{n+1}=b^{n+1},n=0,1,\cdots ,N_t,$$

$$\frac{d{\mu }_{0\mid 1:N_t}}{d{\mu }_0}({\phi }_0)=\frac{d\mu (b^{1:N_t}\mid {\phi }_0)}{d\mu (b^{1:N_t})}.$$

$${\pi }_{0\mid 1:N_t}\propto \pi (b^{1:N_t}\mid {\phi }_0){\pi }_0.$$

$${\hat{\phi }}_{n+1|n}=\mathcal{M}({\phi }_n).$$

$$\frac{d{\mu }_{n+1|n+1}}{d{\hat{\mu }}_{n+1|n}}(\phi )=\exp (-\frac{1}{2}|\mathcal{K}(\phi )-b^{n+1}{|}_{\Gamma }^2).$$

$${\hat{\phi }}_{n+1}^{(j)}=\mathcal{M}({\phi }_n^{(j)}),j=1:J.$$

$${\hat{\Xi }}_n=\frac{1}{J}\sum _{j=1}^J({\hat{\phi }}_{n+1}^{(j)}-{\bar{\phi }}_{n+1})\otimes ({\hat{\phi }}_{n+1}^{(j)}-{\bar{\phi }}_{n+1}),$$

$$\mathfrak{J}(\phi )$$

$$+\frac{1}{2}|\mathcal{K}(\phi )-b^{n+1}{|}_{\Gamma }^2,j=1:J.$$

$${\hat{\varsigma }}_{n+1}=\Theta ({\varsigma }_n)$$

$$\begin{array}{rl}& {\hat{\phi }}_{n+1}=\mathcal{M}({\phi }_n),\\ & {\hat{f}}_{n+1}=\mathcal{G}({\hat{\phi }}_{n+1}).\end{array}$$

$$\mathfrak{J}(\varsigma )=\frac{1}{2}\Vert \varsigma -{\hat{\varsigma }}_{n+1}{\Vert }_{{\hat{\Xi }}_n^{\varsigma }}^2+\frac{1}{2}|\mathcal{O}\varsigma -b^{n+1}{|}_{\Gamma }^2,$$

$$\varsigma$$

$$=(I+{\hat{\Xi }}_n^{\varsigma }{\mathcal{O}}^{\ast }{\Gamma }^{-1}\mathcal{O}{)}^{-1}({\hat{\varsigma }}_{n+1}+{\hat{\Xi }}_n^{\varsigma }{\mathcal{O}}^{\ast }{\Gamma }^{-1}{b}^{n+1}).$$

$$\displaystyle\hat{\Xi}_{n}^{\varsigma}=\left(\begin{array}[]{cc}\hat{\Xi}_{n}^{\varphi\varphi}&\hat{\Xi}_{n}^{\varphi f}\\ \hat{\Xi}_{n}^{f\varphi}&\hat{\Xi}_{n}^{ff}\end{array}\right).$$

Full text

Key words: Level set; Data assimilation; Acoustic source; EnKF

MSC 2010: 35R20, 65R20

1 Introduction

The unknown source problems occur widely in many real applications, e.g., remote sensing, radar detection [24], brain image etc. In this paper, we consider a reconstruction problem of acoustic sources from remote measurements of the acoustic field. For an acoustic source with density $\mathcal{S}(k,x)$ supported on region $D\subset\subset\Omega$, the pressure $u$ of the radiated time-harmonic wave can be modeled [13] by

[TABLE]

where $k=\omega/c_{0}$ is the wave number, $\omega$ is the radial frequency, $c_{0}$ is the speed of sound, $d$ is the spatial dimension and (1.2) is the Sommerfeld radiation condition. Here, we discuss the case of $d=2$ and the source function $\mathcal{S}$ is separable with respect to space and frequency, i.e.,

[TABLE]

In real settings, the source $f$ is usually unknown and needs to be determined from some observed data. And the data receiver devices are generally located on a remote closed surface $\partial\Omega$. We can collect the information of $u$ on $\partial\Omega$ for many frequencies $\omega$ or wave numbers $k$. For $x\in\partial\Omega$, we have the potential representation of acoustic pressure [13]

[TABLE]

where $H_{0}^{(1)}$ is the cylindrical Hankel function. The data are collected on finite points $\vec{x}:=\{x_{1},x_{2},\cdots,x_{N}\}\subset\partial\Omega$ for finite wave numbers $\{k_{1},k_{2},\cdots,k_{M}\}\subset[k_{\min},k_{\max}]$. In Figure 1, we show the measure points by the circles on the boundary of the square. By virtue of (1.4), we can write

[TABLE]

where $\mathcal{H}_{k}:L^{2}(D)\rightarrow\mathbb{C}^{N}$ depends on wave number $k$ and $\eta_{m}$ is the noise. Denote $\mathcal{H}:=[\mathcal{H}_{k}]_{k=k_{1}}^{k_{M}}$ and $b=[b_{1};b_{2};\cdots;b_{M}]$. Further assume that the function $f$ is known a priori to have the form

[TABLE]

here $\mathbb{I}_{D}$ denotes the indicator function of subset $D\subset\mathbb{R}^{2}$, $\{D_{l}\}_{l=1}^{n}$ are subsets of $D$ such that $\bigcup_{l=1}^{n}\bar{D}_{l}=\bar{D}$ and $D_{l}\bigcap D_{j}=\varnothing$ ($l\neq j$), the $\{w_{l}\}_{l=1}^{n}$ are known positive constants. In this setting, the regions $D_{l}$ determine the unknown function and therefore become the primary unknowns. This problem has been discussed in [2, 13], from which we can see that the source $f$ can be determined uniquely according to the multiple frequency data. Many very successful approaches have been proposed to solve the reconstruction problem of the acoustic source [3, 14, 15, 26, 32, 35, 41, 43]. The main techniques in these works include iteration and sampling algorithms. And coping with the ill-posedness is usually addressed through the use of regularization.

For this kind of geometry shape problems, the other very powerful tool is the level set method [29], which was originally introduced by Osher and Sethian [30]. It is devised as a simple and versatile method for computing and analyzing the motion of an interface. In this method, the interface is supposed to move in the normal direction with some speed, and then a partial differential equation, in particular a Hamilton-Jacobi equation, about the level set function is established. Due to the high efficiency, it has developed to be one of the most popular tools for many disciplines, such as image processing, computer graphics, computational geometry. Besides the aspect of high efficiency in computational schemes, the crucial advantage of the level set approach lies in that it allows for topological changes to be detected during the course of algorithms, which is impossible with classical methods based on curve parameterizations. For shape reconstruction inverse problems, level set method has also been widely studied [1, 6, 10, 20, 37].

Recently, statistical perspective has been widely applied in the field of inverse problems [21, 36, 39]. In this framework, the unknown, noise and data are viewed as random variables. Moreover, the core issue for inverse problems in the statistics framework is how to characterize the posterior distribution, e.g., evaluate posterior moments or other posterior expectations. Based on the Bayes’ formula, plenty of algorithms have been proposed to characterize the posterior distribution. The most widely used class of approaches include Markov chain Monte Carlo (MCMC) [36, 40], maximum a posterior estimation (MAP) [36], ensemble Kalman filter (EnKF) [17] etc. These methods has been applied in many inverse problems, e.g., the electrical impedance tomography [12, 22, 23, 34], the Darcy flow [5, 19], fluid mechanics [9, 17], optical flow [4], inverse scattering [7] etc. A natural idea to solve the interface reconstruction problem is to combine the statistical perspective and the level set approach. To the authors’ knowledge, the studies in this field are rare. Since the level set mapping is not continuous, it is difficult to establish the well-posedness of the Bayesian level set posterior distribution according to the theory in [36]. To get the well-posedness, the Gaussian prior assumption on the level set function is imposed in [16] and the almost surely continuous of the forward operator on the level set function is obtained. And subsequently the posterior well-posedness is given [16]. In [11], the hierarchical level set Bayesian algorithm is proposed to solve the kind of problems with unknown prior parameter. In [8], an EnKF coupling level set is considered and applied in EIT problem. One can refer to [27, 28, 31, 38] for some more real applications, e.g., oil reservoir problem, image segmentation problem etc.

2 Level set approach

In this section, we sketch the level set approach for inverse shape reconstruction problems.

According to [10, 11, 16], we can characterize the region $D_{l}$ by the following way

[TABLE]

where $\varphi$ is called the level set function and $-\infty=c_{0}<c_{1}<\cdots<c_{n}=\infty$ are constants. Using this characterization and the assumption (1.6), we can define the level set map $\mathcal{G}:C(\bar{D})\rightarrow L^{2}(D)$

[TABLE]

It is clear that the mapping $\mathcal{G}$ is not an injection since there exist many different $\varphi$ determine the same shape. However, every given function $\varphi$ uniques specifies a function $f$ [10]. We combine (2.2) and (1.5) to give

[TABLE]

This can be rewritten in a compact form

[TABLE]

where $\mathcal{K}=[\mathcal{H}_{k_{m}}\circ\mathcal{G}]_{m=1}^{M}$ and $b=[u(k_{m},x_{j})]$ is an $N\times M$ dimensional column vector.

Introduce a time variable $t$ and define

[TABLE]

The level set function $\varphi$ evolves dynamically given by the convection equation [6, 37]

[TABLE]

which reflects that $\varphi$ varies along with the direction of the negative gradient with velocity $v$. According to [37], we take the variate $\mathfrak{v}$ as

[TABLE]

where $J(f)$ is the Jacobian of $\mathcal{H}f$ at $f$. For an initial guess $\varphi_{0}$, we have the following initial value problem for $\varphi(x,t)$ of Hamilton-Jacobi system

[TABLE]

An explicit discretized version of the above system yields

[TABLE]

with initial guess $\varphi_{0}$.

3 EnKF for the level set algorithm

We now treat problem (2.4) in the frame of statistical inversion, which looks all involved quantities as random variables and explores the posterior distribution on $\varphi$ given $b$. We restate the inverse problem as an identification problem of system state. The system evolves in dynamics (2.10) with unknown initial guess. We assume the initial state $\varphi_{0}$ obeys the distribution $\mu_{0}=N(\bar{\varphi},\Xi_{0})$, where $\bar{\varphi}$ is the prior mean and $\Xi_{0}$ is the prior covariance operator. Observations $b$ are collected at time $t_{n}$ denoted by $b^{n}$. Thereby, we have the following state identification problem

[TABLE]

where $\{\eta^{n+1}\}$ are the i.i.d. Gaussian noise, $\eta^{1}\sim N(0,\Gamma)$ and $\Gamma$ is the noise covariance matrix. We focus on the posterior distribution $\mu_{0\mid 1:N_{t}}$, the distribution of $\varphi_{0}|b^{1:N_{t}}$, $\varphi_{0}$ given $b^{1},b^{2},\cdots,b^{N_{t}}$. Bayes rule gives a characterization of $\mu_{0\mid 1:N_{t}}$ via the ratio of its density with respect to that of the prior

[TABLE]

And therefore, the posterior density satisfies

[TABLE]

From the angle of filter, we are concerned with obtaining the posterior distribution $\mu_{n+1|n+1}$ via $\mu_{n|n}$, the distribution associated with the probability measure on the random variable $\varphi_{n}|b^{n}$. Suppose that we get the present system state $\varphi_{n}$ together with the prior covariance $\Xi_{n}$. The operators $\Xi_{n}$ characterize model uncertainty and are design parameters. With these information as the prior for the next state $\varphi_{n+1}$, we estimate the new state $\varphi_{n+1}$ by blending the observed data $b^{n+1}$ into the dynamics. However, there is, in general, no easily usable closed form expression for the density of $\mu_{n+1|n+1}$. We can break this into a two-step process: prediction step and analysis step. In the prediction step, we use the present state distribution to predict the new state according to the dynamics as

[TABLE]

Using the prediction $\hat{\varphi}_{n+1|n}$ together with the design parameter $\Xi_{n}$, we give the prior distribution $\hat{\mu}_{n+1|n}$ and further obtain the analysis state distribution according to the observational process (3.3)

[TABLE]

When the design parameter $\Xi_{n}$ is fixed, this filter corresponds to the 3DVar process. Naturally, we also can estimate the prior covariance $\Xi_{n}$ via samples. To characterize the posterior distribution, the ensemble Kalman filter (EnKF) approach represents the filtering distribution through an ensemble of particles. The prediction sample can be obtained for each particle. Then the prediction sample mean is taken as the prior mean and the covariance of prediction samples is taken as the prior covariance. In details, suppose that at the $n$-th step, we have samples $\{\varphi_{n}^{(j)}\}_{j=1}^{J}$, $J$ is the number of samples. For each particle, we have a prediction

[TABLE]

Define the empirical covariance $\hat{\Xi}_{n}$ as

[TABLE]

where $\bar{\varphi}_{n+1}=\frac{1}{J}\sum_{j=1}^{J}\hat{\varphi}_{n+1}^{(j)}$ denotes the ensemble mean. The analysis state for each particle is then obtained by minimizing the following functional

[TABLE]

Due to the nonlinearity of the observational operator $\mathcal{K}$ in (3), the numerical implementation of the minimization becomes difficulty. To overcome this point, iterative EnKF approaches are proposed in [17, 18, 42]. Define the state vector $\varsigma=[\varphi,f]^{T}$. We modify the prediction process as

[TABLE]

by

[TABLE]

Accordingly, the observation operator is given by $\mathcal{O}=[0,\mathcal{H}]$. Formally, the minimization function about the state vector $\varsigma$ has the form of

[TABLE]

where $\hat{\Xi}_{n}^{\varsigma}$ is the prior covariance about $\varsigma$. The minimizer can be given by

[TABLE]

Write the prior covariance matrix $\hat{\Xi}_{n}^{\varsigma}$ as

[TABLE]

The component form of (3.14) has

[TABLE]

We list this algorithm here.

Algorithm 1:

• $\bullet$

Initial particles: Let $\{\varphi^{(j)}_{0}\}_{j=1}^{J}$ be the initial samples drawn from the prior distribution $\mathcal{N}(\bar{\varphi},\Xi_{0})$. And obtain the initial samples of $f^{(j)}_{0}$ according to $f^{(j)}_{0}=\mathcal{G}(\varphi^{(j)}_{0})$.

• $\bullet$

Prediction: Evaluate

[TABLE]

Compute $\bar{\varphi}_{n}=\frac{1}{J}\sum_{j=1}^{J}\hat{\varphi}^{(j)}_{n}$, $\bar{f}_{n}=\frac{1}{J}\sum_{j=1}^{J}\hat{f}^{(j)}_{n}$.

• $\bullet$

Analysis: Define $\hat{\Xi}_{n}^{\varphi f}$, $\hat{\Xi}_{n}^{ff}$ by

[TABLE]

Update each particle via (3.18) and (3.19). Compute the mean of the parameter update

[TABLE]

and moreover check the stopping rule.

In other way, the state vector is defined by $\zeta:=[\varphi,\mathfrak{b}]^{T}$ [17]. Accordingly, the observation operator is given as $[0,I]$. In this setting, the EnKF approach is listed in the follows.

Algorithm 2:

1. $\diamond$

Initial particles: Let $\{\varphi^{(j)}_{0}\}_{j=1}^{J}$ be the initial samples drawn from the prior distribution $\mathcal{N}(\bar{\varphi},\Xi_{0})$. 2. $\diamond$

Prediction: Evaluate

[TABLE]

Compute $\bar{\varphi}_{n}=\frac{1}{J}\sum_{j=1}^{J}\hat{\varphi}^{(j)}_{n}$, $\bar{\mathfrak{b}}^{n}=\frac{1}{J}\sum_{j=1}^{J}\hat{\mathfrak{b}}_{(j)}^{n}$. 3. $\diamond$

Analysis: Define $\hat{\Xi}_{n}^{\mathfrak{b}}$, $\hat{\Xi}_{n}^{\mathfrak{b}\mathfrak{b}}$ by

[TABLE]

Update each particle by

[TABLE]

Compute the mean of the parameter update

[TABLE]

and moreover check the stopping rule.

4 Prior distribution

From the above description, we see that the prior distribution $\mu_{0}=N(\bar{\varphi},\Xi_{0})$ plays key role in the statistical inversion. In this paper, we choose the Whittle-Matérn fields as the prior distribution. They are stationary Gaussian random fields with the autocorrelation function

[TABLE]

where $\nu>0$ is the smoothness parameter, $K_{\nu}$ is the modified Bessel function of the second kind of order $\nu$ and $l$ is the length-scale parameter. The integer value of $\nu$ 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 [25, 33, 34]. A method for explicit, and computationally efficient, continuous Markov representations of Gaussian Matérn fields is derived [25]. The method is based on the fact that a Gaussian Matérn field on $\mathbb{R}^{d}$ can be viewed as a solution to the stochastic partial differential equation (SPDE)

[TABLE]

where $W$ is the Gaussian white noise and the constant $\alpha$ is

[TABLE]

Here $\sigma^{2}$ is the variance of the stationary field. The operator $(I-l^{2}\triangle)^{(\nu+1)/2}$ is a pseudo-differential operator defined by its Fourier transform. When $\nu\in\mathbb{Z}$, we can solve (4.2) by the finite element method with suitable boundary condition. The stochastic weak solution of the SPDE (4.2) is found by requiring that

[TABLE]

where $\psi\in L^{2}(D_{1})$. We use the linear element to solve the weak formulation with Dirichlet boundary condition, i.e., $\varphi\mid_{\partial D_{1}}=0$.

5 Numerical test

In the numerical tests, we consider frequencies varying from $k_{\min}=50$ Hz to $k_{\max}=10$ kHz and $c_{0}=343$ms*-1*. These parameter settings are given in [13]. Figure 1 displays the settings of our problem. The data are taken on the circle dots of the boundary $\partial\Omega$. The noise $\eta$ is taken as $\delta*\text{randn}(0,I)$. The prior samples are generated in the disk domain $D_{1}$. We assume that the support $D$ lies in the disk. The disk domain is decomposed into triangular elements and the potential integral is calculated using the Gaussian integral formulation on each triangle. Moreover we solve (4.2) in the same mesh using the linear finite element. To avoid inverse crime as far as possible, we use the finer mesh in the inverse process than in the data generation. The time derivative in the Hamilton-Jacobi system is discretized by the forward Euler scheme. The Euler step size is taken as $0.1$. This sample size $J$ and the iteration step are taken as $1000$ and $100$ respectively. The prior parameter $\nu$ is set to $1$. We display the numerical reconstruction of $f$ for several different cases. In Figures 2, 3 and 4, we give the numerical results for the cases of single domain, two separately domains and Taichi graph respectively. For each case, the data noise level are taken as $0.1$ and $0.01$. It should be noted that we can use (3.24) or (3.25) to generate the estimation of the unknown source $f$ when applying Algorithm 1. We show the numerical comparisons for using different estimation way in Algorithm 1 and Algorithm 2. And it can be seen that whatever approximation way, we almost get similar numerical result for all examples. In addition, it is obvious that the numerical effect is better for smaller noise level.