Bayesian approach for inverse obstacle scattering with Poisson data

TL;DR

This paper develops a Bayesian inversion method for reconstructing acoustic obstacles from Poisson data, incorporating hybrid priors including Gaussian and total variation to improve the reconstruction quality.

Contribution

It introduces a Bayesian framework with hybrid priors for inverse obstacle scattering with Poisson data, and discusses the well-posedness of the posterior distribution.

Findings

Numerical examples demonstrate the effectiveness of the proposed algorithm.

Hybrid priors improve reconstruction accuracy.

The approach handles stochastic Poisson data effectively.

Abstract

We consider an acoustic obstacle reconstruction problem with Poisson data. Due to the stochastic nature of the data, we tackle this problem in the framework of Bayesian inversion. The unknown obstacle is parameterized in its angular form. The prior for the parameterized unknown plays key role in the Bayes reconstruction algorithm. The most popular used prior is the Gaussian. Under the Gaussian prior assumption, we further suppose that the unknown satisfies the total variation prior. With the hybrid prior, the well-posedness of the posterior distribution is discussed. The numerical examples verify the effectiveness of the proposed algorithm.

Equations

$$△u+k^{2}u=0\text{in}{\mathbb{R}}^2\backslash D,$$

$$u{\mid }_{\partial D}=0,$$

$$\underset{r\to \infty }{\lim }\sqrt{r}(\frac{\partial u^{\mathrm{s}}}{\partial r}-ik{u}^{\mathrm{s}})=0,r=|x|.$$

$$u^{\mathrm{s}}(x)=\frac{{\mathrm{e}}^{ik|x|}}{|x{|}^{\frac{1}{2}}}\left\{u_{\infty }(\hat{x})+\mathcal{O}\left(\frac{1}{|x|}\right)\right\}\text{as}|x|\to \infty ,$$

$$F(D)=\tau |u_{\infty }{|}^2,$$

$$\partial D=q(t)[\cos (t),\sin (t){]}^T,t\in [0,2\pi ).$$

$$G(q)=\tau |u_{\infty }{|}^2.$$

$$(S\phi )(x)=2{\int }_{\partial D}\Phi (x,\tilde{x})\phi (\tilde{x})ds(\tilde{x}),x\in \partial D,$$

$$(K\phi )(x)=2{\int }_{\partial D}\frac{\partial \Phi (x,\tilde{x})}{\partial \nu (\tilde{x})}\phi (\tilde{x})ds(\tilde{x}),x\in \partial D,$$

$$u^{\mathrm{s}}(x)={\int }_{\partial D}\left\{\frac{\partial \Phi (x,\tilde{x})}{\partial \nu (\tilde{x})}-i\eta \Phi (x,\tilde{x})\right\}\phi (\tilde{x})ds(\tilde{x}),x\in {\mathbb{R}}^2\backslash \overline{D},$$

$$(I+K-i\eta S)\phi =-2u^{\mathrm{i}}\text{on}\partial D.$$

$$u_{\infty }(\hat{x},d)=\frac{e^{-i\frac{\pi }{4}}}{\sqrt{8\pi k}}{\int }_{\partial D}(k\nu (\tilde{x})\cdot \hat{x}+\eta )e^{-ik\hat{x}\cdot \tilde{x}}\phi (\tilde{x})ds(\tilde{x}).$$

$$\frac{d{\mu }^y}{d{\mu }_{\mathrm{pr}}}=\pi (y|q),$$

$$\pi (y|q)=\prod _{j=1}^m\frac{{\lambda }_j^{y_j}{e}^{-{\lambda }_j}}{y_j!}.$$

$$\pi (y|q)\propto \exp (-\Lambda (q,y)),$$

$$\Lambda (q,y)=\sum _{j=1}^m(G(q){)}_j-y_j\log (G(q){)}_j.$$

$$q(t)=\exp (r(t)),$$

$$q(t)=\frac{a}{2}(\mathrm{erf}(r(t))+b),$$

$$\mathrm{erf}(r(t))=\frac{2}{\sqrt{\pi }}{\int }_0^{r(t)}\exp (-t^2)dt.$$

$${\mu }_{\mathrm{pr}}={\mu }_1=N(0,C_1),$$

$$r^{\prime \prime }(t)\sim N(0,A^{-s}),s>0,$$

$$D(A):=\left\{v\in H^2[0,2\pi ]:{\int }_0^{2\pi }v(t)dt=0\right\}.$$

$$C_2(s,t)=\exp (-\frac{2{\sin }^2(\frac{t-s}{2})}{l^2}),$$

$$|r(t)-r(s){|}^2=(\cos t-\cos s{)}^2+(\sin t-\sin s{)}^2=4{\sin }^2\frac{t-s}{2}.$$

$$C_2(r(t),r(s))=\exp (-\frac{|r(t)-r(s){|}^2}{2l^2}).$$

$$\frac{d{\mu }_{\mathrm{TV}}}{d{\mu }_2}(r)\propto \exp (-\Vert r{\Vert }_{\mathrm{TV}}):=\exp (-R(r)),$$

$$\Vert r{\Vert }_{\mathrm{TV}}=\zeta {\int }_0^{2\pi }|\frac{dr}{dt}|dt$$

$$\frac{d{\mu }^y}{d{\mu }_2}\propto \exp (-\Lambda (r,y)-R(r)):=\exp (-\Psi (r,y)).$$

$$d_{\mathrm{H}}({\mu }^y,{\mu }^{\tilde{y}})={\left(\frac{1}{2}\int (\sqrt{\frac{d{\mu }^y}{d{\mu }_2}}-\sqrt{\frac{d{\mu }^{\tilde{y}}}{d{\mu }_2}}{)}^{2}d{\mu }_2\right)}^{\frac{1}{2}}.$$

$${\int }_X\exp (\alpha \Vert r{\Vert }_X^2)d\mu <\infty .$$

$$|G(r)|⩽M(1+\Vert r{\Vert }_X).$$

Full text

Key words: Bayesian inversion; Acoustic scattering; Poisson data; Total variation prior; MCMC

MSC 2010: 35R20, 65R20

1 Introduction

Inverse obstacle scattering problems occur in many real-world applications such as radar, sonar and non-destructive testing [2]. One main object of this kind of problems is to find the shape and location of the unknown obstacle, which usually can not be observed directly. Due to the widely applications, these prolems have attracted many researchers’ attention. Especially, a large number of papers put emphasis on designing some effect numerical algorithms to recover the unknown information according to some related observed data [1, 2]. The observed data are contaminated by the noise inevitably. The Gaussian process has many theoretical and practical advantages since it can mimic the effect of many random processes. And therefore a lot of noise models can be characterized by the Gaussians. However, in many applications, the response variable of interest is a count, that is, takes on nonnegative integer values. For count data, it is well modeled by a Poisson process. The Poisson process is usually used to describe the occurrences of events and the time points at which the events occur in a given time interval, such as the occurrence of natural disasters and the arrival times of customers at a service center. In the obstacle scattering problem, the counts of photons at a remote closed surface are detected. As seen in [3, 9, 10], the counts of photons have interacted with the unknown object and can be well modeled by a Poisson process. In view of this, we consider to use the Poisson data to recover the unknown shape.

The stochastic nature of the data undoubtedly results in the uncertainty of the reconstruction. Therefore, besides the reconstruction of the unknown object itself, its uncertainty requires the extra concerns. Probabilistic thinking provides a natural way to quantify the uncertainty. Recently, Bayesian inference method has become a popular tool for this purpose [7]. In Bayesian approach, the main idea is to translate the prior knowledge on the noise and on the unknowns to prior probability laws and then use the forward model to transmit the prior information to the posterior distribution from which we can deduce any information about the unknowns. A rigorous Bayesian framework is developed for the inverse problems with the infinite dimensional unknowns [8]. From the existed works, we can see that the prior for the unknowns is crucial in Bayesian approach and even determines whether or not it works. A difficult subject is how to code the experience prior information in some distribution formula. It is well-known that, subject to certain mild conditions, the sum of a set of random variables, has a distribution that becomes increasingly Gaussian as the number of terms in the sum increases. This point allows us to use Gaussians to model the prior for many real physical parameters. In this paper, we characterize the unknown by two kinds of prior knowledges: the Gaussian and the hybrid priors. Upon these priors, we discuss the well-posedness of the posterior distribution with the Poisson data.

The remainder of this paper is organized in the following: In Section 2, we describe the exterior acoustic scattering problem. In Section 3, the Bayesian approach and the well-posedness are discussed. In the final, some numerical examples are given.

2 Acoustic obstacle scattering with Poisson data

For a perfect cylindrical conductor with smooth cross section $D\subset\mathbb{R}^{2}$, the scattering of polarized, transverse magnetic time harmonic electromagnetic waves is described by the Helmholtz equation

[TABLE]

where $k>0$ is the wave number, $u=u^{\rm s}+u^{\rm i}$ the total field, $u^{\rm s}$ the scattered field and $u^{\rm i}$ the incident wave. We let $u^{\rm i}=\exp(ikx\cdot d)$, where $d$ is the incident direction. The following boundary conditions are imposed

[TABLE]

The condition (2.2) corresponds to the sound soft boundary condtion and (2.3) is the so-called Sommerfeld radiation condition. The Sommerfeld radiation condition characterizes the outgoing wave and enables us to have the following asymptotic behaviour for the scattered wave $u^{\rm s}$:

[TABLE]

where $\hat{x}=\frac{x}{|x|}$ and $u_{\infty}(\hat{x})$ is the far field pattern. The direct scattering problem is to find the scattered field $u^{\rm s}$ for given incident field $u^{\rm i}$ and obstacle $D$. It is known that there exists a unique solution $u\in H^{1}(\mathbb{R}^{2}\backslash D)$ if $\partial D$ is Lipschitz continuous [2].

We focus on the reconstruction of the unknown domain $D$ by some observed data: the photon counts of the scattered electromagnetic field far away from the obstacle. Since at large distances the photon density is approximately proportional to $|u_{\infty}|^{2}$ [9], the obstacle reconstruction problem can be described by the operator equation

[TABLE]

where $\tau$ is a constant corresponding the noise level. We suppose that the domain $D$ is starlike, i.e., the boundary $\partial D$ is parameterized by

[TABLE]

In this setting, the function $q$ becomes the unknown and therefore the equation (2.5) is equivalent to

[TABLE]

The operator equation can be given by the single-double layer potential theories [2]. Let $\partial D$ be of class $C^{2}$. Define the single layer potential operator $S$

[TABLE]

and the double layer potential operator $K$

[TABLE]

where $\Phi(x,\tilde{x})=\frac{i}{4}H_{0}^{1}(|x-\tilde{x}|)$ is the fundamental solution, $H_{0}^{1}(\cdot)$ is the Hankel function of the first kind of order zero and $\nu(\tilde{x})$ is normal direction. It can be seen in [2] that $S$ and $K$ are bounded from $C^{0,\alpha}(\partial D)$ into $C^{1,\alpha}(\partial D)$. With the single-double layer potentials, the scattered field can be written as

[TABLE]

where $\eta$ is the real coupling parameter and $\varphi$ the unknown density function. Then the main goal for the direct scattering problem is to determine the unknown density function $\varphi$ according to the sound soft boundary condition by

[TABLE]

By (2.11), the far field pattern further can be written as

[TABLE]

Equations (2.11), (2.12) together with (2.6) give the forward operator $G$.

3 Bayesian inversion

Denote $X=C^{2,\alpha}[0,2\pi)$ and $Y=\mathbb{R}^{m}$. We consider the observation $y\in Y$ is a Poisson point process [9] in $[0,2\pi)$ with parameters $\lambda_{j}=(G(q))_{j}=(\tau|u_{\infty}|^{2})_{j},\,\,j=1,2,\cdots,m$. The goal of statistical inversion is to explore the posterior distribution $\mu^{y}(q)$, i.e., the distribution of $q$ conditional on the value of $y$. Bayes’ formula shows that

[TABLE]

where $\mu_{\rm{pr}}$ is the prior distribution for the unknown $q$ and $\pi(y|q)$ the likelihood function.

The Poisson noise yields the likelihood function

[TABLE]

Following [7, 8], the likelihood function can be written in the form of

[TABLE]

where

[TABLE]

The prior is given before acquiring the data. Since the $q$ denotes the length, an additional restriction should be imposed: the unknown must be positive. For this purpose, we consider the following two ways [1, 10]

[TABLE]

where $a>0$ and $b>1$ are constants and $\rm{erf}(\cdot)$ is the error function

[TABLE]

These two ways guarantee that $q$ is strictly positive. With these parameterizations, we now infer the new unknown parameter $r$. We still use $\mu_{\rm pr}$ to denote the prior distribution for $r$. Without any confusion, the forward operator is still denoted by $G$. The prior needs to be chosen carefully and it is not unique. We give different priors for (3.5) and (3.6).

For (3.5), the prior $\mu_{\rm pr}$ is given by the Gaussian

[TABLE]

where $C_{1}$ is the covariance function and should be chosen carefully. In [1], a convenient prior to use is

[TABLE]

where $A=\frac{d^{2}}{dt^{2}}$ with the definition domain

[TABLE]

In this prior assumption, it can be verified $r\in X$ almost surely with $0<\alpha<\min\{1,s-1/2\}$, i.e., $\mu_{1}(X)=1$ [1, 7]. In this prior assumption, the well-posedness of the posterior distribution for the Gaussian noise case has been discussed in [1].

For (3.6), the prior distribution is chosen as $\mu_{\rm pr}=N(0,C_{2})$ denoted by $\mu_{2}$. To guarantee the period, we use the covariance function [5, 6]

[TABLE]

where $l$ defines the characteristic length-scale. The covariance example can be viewed as mapping the one-dimensional input variable $t$ to the two-dimensional $r(t)=(\cos t,\sin t)$. The distance between $r(t)$ and $r(s)$ is given

[TABLE]

The covariance kernel (3.10) can be obtained by taking the squared exponential kernel in $r$ space

[TABLE]

We apply a hybrid prior proposed in [10]. On the basis of the Gaussian prior $\mu_{2}$, we further assume that the unknown $r$ satisfies the total variation prior as in [10]

[TABLE]

where $\|\cdot\|_{\rm TV}$ is the TV seminorm,

[TABLE]

with a positive constant $\zeta$. It follows immediately that the posterior distribution has this form

[TABLE]

Next, we discuss the well-posedness of the posterior distribution with the hybrid prior. This well-posedness is characterized in the sense of Hellinger metric, which is defined by

[TABLE]

Here, we need to use the following Fernique theorem to some functionals are bounded for Gaussian measures.

Theorem 3.1** (Fernique).**

If $\mu=N(0,C)$ is a Gaussian measure on Banach space $X$, so that $\mu(X)=1$, then there exists $\alpha>0$ such that

[TABLE]

Lemma 3.2*.*

For fixed $\hat{x},d$, there exists $M>0$ such that

[TABLE]

for all $r\in X$.

Proof*.*

According to (2.6) and (2.12), we have

[TABLE]

where $\tilde{\varphi}(t)=\varphi(x)\mid_{x=x(t)}$. Hence it follows that

[TABLE]

By (3.6), we get

[TABLE]

This completes the proof.

Lemma 3.3*.*

For every $s>0$, there exists constants $I_{2}(s)\in\mathbb{R}$ and $I_{1}(s)>0$ such that, for all $r\in X$, and $y\in Y$ with $|y|<s$,

[TABLE]

Proof*.*

For convenient, denote

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the inner product in $\mathbb{R}^{m}$. It is obvious from (3.6) that

[TABLE]

Therefore, the boundedness of the forward operator $G$ implies that

[TABLE]

for two constant vectors $\underline{\lambda}$, $\overline{\lambda}$. And it follows immediately that $|\log\lambda|\leq l_{\rm max}$ for a positive constant $l_{\rm max}$. Furthermore, we obtain for $|y|<s$ by the Cauchy-Schwarz inequality

[TABLE]

Remark 1*.*

It should be noted that the lower bound $\underline{\lambda}$ in (3.16) probably has [math] component. In this case, we can make a shift for the forward operator by

[TABLE]

where $e$ is a positive constant vector. Then we discuss the new problem of the form

[TABLE]

To avoid the addition statement, we still use $G$ as the forward operator.

Theorem 3.4**.**

Assume that $\mu_{2}$ is the Gaussian measure with covariance function (3.10). For $y,\tilde{y}$ with $\max\{|y|,|\tilde{y}|\}<\varpi$, there exists $M=M(\varpi)>0$ such that

[TABLE]

Proof*.*

First, it should be pointed out that the posterior measure $\mu^{y}$ is absolute continuous about the prior measure $\mu_{2}$. This will guarantee the existence of the Radon-Nikodym derivative $\frac{d\mu^{y}}{d\mu_{2}}$. In fact, for any measurable subset $X_{1}\subset X$, it holds that

[TABLE]

Therefore $\mu_{2}(X_{1})=0$ implies that $\mu^{y}(X_{1})=0$. Define

[TABLE]

It is obvious that

[TABLE]

By Lemma 3.3, we have

[TABLE]

Since $\mu_{2}$ is the Gaussian measure, the unit ball has positive measure. This shows that $A(y)$ is strictly positive. This together with this absolute continuous shows that the posterior distribution $\mu^{y}$ is well-defined. The Hellinger metric holds that

[TABLE]

For $\Sigma_{1}$, since $A(y)$ is bounded, it suffices to estimate the integral part. According to Lemma 3.2 and Fernique therorem, we have

[TABLE]

Since $A(y)$ and $A(\tilde{y})$ are bounded from below, it follows that

[TABLE]

The estimation (3.17) can be obtained from (3.18), (3.19) and (3.20).

Remark 2*.*

The requirement $\mu_{2}(X)=1$ in Fernique theorm is guaranteed naturally since it has the exponential covariance function (3.10). This case for $q$ with the form (3.5) is similar to (3.6). Actually, we just need to take $\zeta=0$ and verify the condition $\mu_{1}(X)=1$.

4 Numerical test

In this section, we demonstrate some numerical examples to show the effectiveness of the Bayesian method for the inverse obstacle scattering problem with the Poisson data. We are ready to use the MCMC method to generate samples to explore the posterior distribution $\mu^{y}$ in (3.1).

It is clear that the operator $A$ in (3.9) has the eigen-system $\{n^{2};\frac{\cos(nt)}{\sqrt{\pi}},\frac{\sin(nt)}{\sqrt{\pi}}\}$. With this, the prior for (3.5) can be generated by the Karhunen-Loève expansion. The truncated expansion is used to generate the prior samples. The corresponding derivatives can be computed accurately. For (3.6), the samples for $r$ with the covariance (3.10) are period. We use the fast Fourier transform (FFT) to implement the numerical computations for the derivatives of the samples. The posterior samples for the two ways are both generated by the Metropolis-Hastings algorithm.

We just list the hybrid prior case in the following. For using (3.5) and the corresponding prior, the process is similar.

• $\bullet$

Initialize $r^{0}$ and further generate $q^{0}$ in (3.6).

• $\bullet$

For iteration $i=1,2,\cdots$ do

• $\cdot$

Propose $\hat{r}=\sqrt{1-\beta^{2}}r^{i-1}+\beta\xi,\,\,\xi\sim N(0,C_{2})$ and further generate $\hat{q}$.

• $\cdot$

Acceptance probability:

[TABLE]

• $\cdot$

Accept or reject the proposal in the following way

[TABLE]

endfor

Consider the following obstacles

[TABLE]

where $t\in[0,2\pi)$.

In Fig. 1, we plot the Poisson data for $\tau=1000$. Using these data, we reconstruct the obstacles displayed in Fig. 2 and 3. The posterior sample means are used as the approximations. In Fig. 2, we plot the $95\%$ confidence interval for the peanut shape.