Bayesian Method for Shape Reconstruction in the Inverse Interior Scattering Problem

The inverse scattering problem of an interior cavity with three different boundary conditions is considered. Bayesian method is used to reconstruct the shape of the cavity from scattered fields incited by point source(s) and measured on a closed curve inside the cavity. We prove the well-posedness in Bayesian perspective and present numerical examples to show the viability of the method.


Introduction
Recently, the inverse scattering problem of an interior cavity attracts many researchers' attention due to its useful applications in industry.For example, in nondestructive testing the sources and receivers are sometimes placed inside the cavity to test the structural integrity [1].Different from the typical inverse obstacle problem which is an exterior boundary value problem where the wave incidence and measurements are taken outside the obstacle, our problem is to recover the shape of the cavity from point source(s) and measurements on a closed curve inside the cavity.To be precise, we denote the cavity by a bounded simply connected domain  ⊂ R 2 with Lipschitz boundary .Our aim is to reconstruct the boundary  from measurements (scattered waves) taken on a closed curve  inside the cavity  (see Figure 1).The point source(s) (incident waves) are also located on curve .Some numerical methods have been proposed for this kind of inverse problem.The first paper related to this interior problem with Dirichlet boundary condition is Qin and Colton [2].They proved the uniqueness and used a modification of linear sampling method to reconstruct the shape of the cavity.Then they developed this method for impedance boundary condition [3] and Hu et al. developed this method for mixed boundary condition [4].In 2011, Qin and Cakoni [5] proposed a nonlinear integral equation method and Zeng et al. [6] developed the linear sampling method for inverse interior electromagnetic scattering problem.A decomposition method has been presented by Zeng et al. in [7].Recently Liu [8] proposed a factorization method for this inverse interior cavity problem.
Bayesian theory is the central part of statistical inversion theory, so sometimes we just call the statistical inversion method Bayesian method.Different from the traditional approach which produces single estimate of the unknowns, Bayesian method produces a distribution which describes the behaviour of solution based on the prior information and the random measurement.All variables in the model are viewed as random variables and the randomness that is the degree of information about these variables is coded in probability distributions.Bayesian approach plays same role as regularization methods when dealing with the ill-posed inverse problems.Every prior distribution can be replaced by an appropriately chosen penalty and it tells us the information hidden in the penalty.Another advantage of Bayesian method is that it is a good interpretation of mathematical models, well understood, and generally accepted [9].From numerical perspective, Bayesian approach is easy to encode when you already have the numerical method for direct problem.The process of solving direct problem can be viewed as a black box, because we are just concerned with the input and output.Although the computational amount is large for we should sample many points to describe the posterior distribution, Bayesian method is already used in many application areas with the development of computational capabilities.For some classical books and papers about Bayesian method we refer to [9][10][11][12][13][14].
Due to the advantages stated above, we choose Bayesian method to solve the inverse interior cavity problem.The study of using Bayesian approach for shape reconstruction in inverse scattering problem is a fairly new topic.It brings the randomness into the deterministic problem and provides us with a new perspective to view our problem.In the traditional method, the reconstructed geometries are different when the input of random noise varies and we do not know which one is more reliable.But in the Bayesian method we can directly obtain the distribution of our reconstructed result and it will not change when the distribution of noise is decided.From numerical simulation of the distribution we know which reconstructed geometries are more reliable (with high probability).Then we can use the statistic characteristics of this distribution, for example, the mean value, to estimate the result we wanted.Comparing to the sampling method, Bayesian method does not need to measure scattered fields corresponding to many point sources, respectively.Only measuring the scattered field once is enough.Of course, there are disadvantages of Bayesian method; for example, the computational efficiency is not high.But it gives us a new way to understand and cope with inverse interior cavity problem and it can be also extended to other inverse scattering problems, such as inverse open cavity problem and inverse rough surface problem.In our paper, we adopt the framework in [14] for an interior inverse scattering problem with three different boundary conditions.We prove the wellposedness of Bayesian method for this problem and present some numerical results.
The outline of this paper is as follows.In Section 2 we describe the problem and give the Bayesian formulation of our problem.In Section 3, we prove two important properties of the direct solution operator and get the well-posedness of Bayesian method.In Section 4, some numerical examples are presented to show the effectiveness of our method.

Bayesian Formulation of the Problem
In this paper we consider a TM polarized time harmonic electric dipole located inside an infinite cylinder.Let the cross section of the cylinder  ⊂ R 2 be a simply connected domain with  2, boundary .We assume the point sources and the observational points all locate on curve  inside the domain .Then the above scattering problem reduces to finding the scattered field   which satisfies where  > 0 is the wave number and   is the incident wave given by Here,  = √ − Here,  > 0 and ] is the unit outward normal to .
The direct problem is solving the above equation to get scattered field   for given .The inverse problem is recovering the shape  for given measurements   .To ensure the uniqueness of interior scattering problem, we assume that  2 is not an eigenvalue for −Δ corresponding to boundary condition (2) of .For simplification, we assume  is a starlike domain with respect to the origin.Thus, we can represent the boundary  by where 0 <  <  0 .We set  = ln just for the convenience of the proof in the next section.Then our model can be written as where G := (  ( 1 ), . . .,   (  )) ( is the number of measurements, i.e., y ∈ R  ) is a finite dimensional observational operator corresponding to (1) and vector y is the observational data with noise .We describe our prior information about  ∈ , in terms of a probability measure  0 .Here,  is some function space to be chosen (see Section 3.1).And we use Bayes' formula to calculate the posterior probability measure   for  given y.Let  0 and   denote the probability density functions of measures  0 and   .By the knowledge in statistics the noise is often a mean zero random variable and we know little about  as prior information in the interior inverse scattering problem.So we assume the noise  ∼ N(0, Γ) with covariance matrix Γ and the prior  0 ∼ N( 0 ,  0 ) with mean function  0 and bounded covariance operator  0 .Of course, the prior measure  0 may be various from the different prior information of the cavity geometry  and we will present a numerical example in Section 4.
According to Bayes' formula, we can get where | ⋅ | 2 Γ := (⋅, Γ −1 ⋅) and  ∝  means  is proportional to .Then the Radon-Nikodym derivative related to the prior measure and posterior measure is In this paper, we use MCMC (Markov Chain Monte Carlo) method to describe the posterior distribution according to   () and get the mean of .For the details of this method we refer to [12,15].

Well-Posedness of Bayesian Method
In this section, our goal is to enable application of the framework in [14] for our inverse interior scattering problem.We divide this section into two subsections.In Section 3.1 we give the well-posedness of our Bayesian method and in Section 3.2 we prove that our observational operator G satisfies the assumption in the well-posedness theorem.

Well-Posedness Theorem.
According to Section 4 (common structure) of [14], we should choose a space  and a prior measure  0 such that  0 () = 1 and Assumption 2.7 in [14] holds on .Then we can get the well-posedness results of Bayesian method.So first we give the specific definition of prior measure on .Assume  = − 2 / 2 with the definition domain: Because in the representation of starlike domain () is 2 periodic, here and in the remaining part of the paper we use the square bracket in the definition of the domain [0, 2] to signify the periodicity.So  2 [0, 2] is the usual Sobolev space with periodic condition (0) = (2); that is, Now we assume a prior Gaussian measure on  such that   () ∼ N(0,  −1 ).Through simple calculation, we can get that the eigenvalues of  −1 are 1/ 2 and the corresponding eigenfunctions are  2 = cos()/√ and  2−1 = sin()/√.According to Karhunen-Loève expansion, we have where   and   are i.i.d.(independent and identically distributed) sequences with  1 ∼ N(0, 1) and  1 ∼ N(0, 1).Clearly, this is in fact the Fourier expansion.Integrating   () we can get ().() is also a Gaussian measure because the integral operator is linear and continuous [16].Using Lemma 6.25 in [14], we know that   () is almost surely in  0, for any  < 1/2.In order to get (), we integrate   ().The function () is not unique when its second derivative   () is given.
Here, we just use periodic form of () in ( 13) the same as in [10].More specifically, we define where  0 ∼ N(,  2 ) is a Gaussian random variable.Then we can integrate the Fourier expansion of   () term by term to obtain Here, we define the Gaussian measure on  based on its second derivative   .There is also an alternative way to define the Gaussian measure on ; see [10].Now we define  = where Then we obtain the following theorem.
Proof.From the analysis above, we know   () is almost surely in  0, for any From Theorem 4.1, Theorem 4.2, Theorem 6.31, and Lemma 2.8 in [14], in order to get the well-posedness theorem, we need to prove that the observational operator G satisfies the following assumption.
Assumption 2 (Assumption 2.7 in [14]).(i) For every  > 0, there is an  = () ∈ R such that, for all  ∈ ,     G ()    Γ ⩽ exp ( (ii) For every  > 0, there is a  = () > 0 such that, for all We will prove that G() satisfies Assumption 2 in the next subsection.Now we give the main theorem in our paper.
(i) The posterior measure   is absolutely continuous with respect to  0 and has Radon-Nikodym derivative given by (ii) The posterior measure   is a well-defined probability measure on  2 [0, 2].(iii) The posterior measure   is Lipschitz in the data y, with respect to the Hellinger distance; that is, there exists a constant () > 0 such that for all y, y  with max{‖y‖ Γ , ‖y  ‖ Γ < }.Here, the Hellinger distance is defined as provided  1 and  2 are both absolutely continuous with  0 .
The proof of Theorem 3 is just a result of application of Theorem 4.1, Theorem 4.2, Theorem 6.31, and Lemma 2.8 in [14].In fact, from Theorem 4.2 in [14] we know the mean of   is also continuous in y.So if we use the mean of posterior distribution to approximate the solution () of inverse interior scattering problem, it is continuous with the observational data y.

Properties of Observational Operator G.
In this subsection, our task is to prove Assumption 2 when G() is defined by integral equation method in [17,18].We use the classical layer potential approach to formulate integral equations for direct interior scattering problem with three different boundary conditions.First, let us introduce the single-and double-layer operators  and , given by and the normal derivative operator   , given by where Φ(, ) is the fundamental solution of Helmholtz equation and ]() is the unit outward normal.
Then we look for a solution   in the form of for the Dirichlet boundary condition, or for the Neumann boundary condition or impedance boundary condition.From the jump relations, we see that   is the solution of scattering problem (1)-( 2), provided the density  is a solution of the following integral equation: for the Dirichlet boundary condition, or for the Neumann boundary condition, or for the impedance boundary condition.
From [19] we known that the operators  and  are compact operators on  2, (), and   is a compact operator on  1, ().If  2 is not an eigenvalue of −Δ for all possible domain  when  changes, we can get the injectivity of  − ,  +   , and  +   + .For details we refer to [17].In practice, the cavity is always bounded, so () is also bounded.We can choose positive constants  min small enough and  max big enough such that  min <  <  max .Then by Riesz-Fredholm theory [17], from the injectivity of the operators  − ,  +   , and  +   + , we know that they are bijective, their inverse operators are bounded on domain , and the bounds of the operators are only related to the constants  min and  max .From the analytic property of   on , we obtain () ∈  2, () for Dirichlet boundary condition and () ∈  1, () for Neumann or impedance boundary condition.S maps  1, () continuously into  2, ().As a result,   () ∈  2, () and the trace   ()|  ∈  2, () for all the three boundary conditions [19].

Remark 4.
Here, for the impedance boundary condition, when  > 0 is a real number,  2 is not an eigenvalue for −Δ because of the uniqueness of interior Helmholtz problem with impedance boundary condition.So the case of impedance boundary condition does not suffer from the existence of eigenfrequencies [5,20].
For the Dirichlet boundary condition, we can find a  small enough to avoid all Dirichlet eigenvalues.More specifically, for  <  max , the ball   max contains all possible domain .From the Rayleigh-Faber-Krahn inequality the smallest Dirichlet eigenvalue  1 (>0) of −Δ for domain  is no less than the smallest eigenvalue   1 (>0) of   max .Then when  2 <   1 , it is not an eigenvalue of −Δ for all possible domain .This method of avoiding eigenvalues is often used in inverse interior scattering problem [5].
For the Neumann boundary condition, we cannot give a similar method to avoid all the eigenvalues as Dirichlet boundary condition.The distribution of Neumann eigenvalues of interior scattering problem is still under investigation.So for Neumann boundary condition we just assume that  2 is not an eigenvalue of −Δ for all possible domain .In fact, if the cavity is filled with some special material such that Im > 0, then the interior Dirichlet and Neumann problems have at most one solution [17].We will consider this situation in the future investigation.(33)

Now we show that G(𝑞) satisfies the property (i) in
Now our goal is to bound the two factors in right hand of the above inequality.
That is, G() satisfies property (ii) in Assumption 2.
In order to get this Lipschitz continuity, we use the definition of domain derivative in [21].There are also other ways to define and get the domain derivative, such as [22][23][24], but the results are similar.Assume  is a bounded open domain in R 2 with  2 smooth boundary .() is a function defined on  satisfying where  and  are some partial differential operators.Given  a regular (in our application regular means  2 smooth) vector field defined in R 2 , one denotes For  small enough, ( + )() is an open set close to  with regular boundary.( + )() can be viewed as a variation of .Then we define () = (( + )()) and we call the local derivative (()/)(0) the domain derivative with .Let   = (()/)(0) be the local derivative of  in a direction ; that is, where  is a regular vector field.
In our problem, the smoothness of  and  and the regularity of  satisfy the hypotheses (3.1) to (3.9) and (3.12) in [21].The operator  = Δ +  2 ,  = B,  = 0, and  = −B  .Clearly,  and  are linear bounded operators, so the derivatives are themselves; that is, Then we obtain the derivative   in our problem satisfies Here, because of the boundary condition (2), the tangential part of ∇(B(  +   )) is 0; that is, In our problem,  is a starlike domain.The variation of  is in fact the variation of .So we can denote  = exp()(cos, sin)q, where q is the direction of variation of .Then the boundary condition (44) changes to where   = (cos, sin).Now we show that G() satisfies property (ii) in Assumption 2. Theorem 6.For every  > 0, there is a  = () > 0 such that, for all  1 ,  2 ∈  with

Numerical Examples
We have already demonstrated the well-posedness of the posterior distribution.Now we give some numerical examples to show the effectiveness of our method.In order to describe the posterior distribution, we consider to adopt MCMC method to generate samples distributed according to (7).Then the average of these samples can be used to approximate the expectations with respect to the posterior distribution and hence to make predictions about the shape parameters [12].Recently, a DRAM (Delayed Rejection Adaptive MCMC) algorithm is proposed to improve the efficiency of the standard MH (Metropolis-Hasting) algorithm.DRAM algorithm combines two powerful ideas in MCMC: adaptive Metropolis samplers and delayed rejection [15].So in this paper, we use DRAM algorithm to generate samples.We consider a peanut cavity which is the radial function parameterized by () = 0.4 √ 4cos 2  + sin 2 .The interior curve  is a circle defined by  = {() | () = 0.3(cos, sin),  ∈ [0, 2]}.We take the source points  1 = (0, 0.3) and  2 = (0, −0.3) on .The wave number  = 1.In order to avoid the inverse crime we solve the direct problem by MFS (method of fundamental solution) to obtain the synthetic data (observational data G()) on curve .Three different boundary conditions are considered, respectively, in Examples 1, 2, and 3.All the numerical experiments followed were performed using MATLAB software.
Example 1 (the sound-soft boundary condition).We truncate 6 terms of the series in (13) and assume the corresponding 13 coefficients obey Gaussian distribution N(0, 1).We take the measurement error  ∼ N(0,  2 ),  = 0.005 (the relative error is about 0.3%), and generate 100000 samples of the posterior distribution.The posterior distribution histograms of coefficients  0 and  1 are presented in Figure 2 and the Markov Chains of them are presented in Figure 3. From Figure 2 we can see the general shape of the posterior distributions of  0 and  1 .We omit the histograms for the other coefficients due to the space constraints.The trace plots of the Markov Chains in Figure 3 look like "fuzzy worms, " so as a rule of thumb the step size in our algorithm is appropriate [12].
Because the beginning draws of the chain represent poorly the distribution to be explored [12], we remove the first 20000 draws from the chain and take the mean value of the other draws.The reconstruction of the cavity is shown in Figure 4.The reconstruction of the radial function () and the scattered field   are displayed in Figure 5, where the -axis represents the angle value varying from 0 to 2.In Figure 4 and left one in Figure 5 the red line is the exact cavity shape and the blue line is the reconstruction.For the right figure of Figure 5 the solid line is the exact value and the "+" line is the reconstruction value.We use the same notations in the following experiments without extra illustrations.From these three reconstruction figures we can see our method is effective.
In the next experiment, we consider the effect of different levels of noise.In Figures 6 and 7, the noise  ∼ N(0,  2 ),  = 0.01 (the relative error is about 0.6%).In Figures 8 and 9, the noise  ∼ N(0,  2 ),  = 0.05 (the relative error is about 3.3%).In the right one of Figure 9, the green line denotes the measurement data y (the exact scattered field added with noise).From the reconstruction of scattered field, we know that, in the interior cavity scattering problems, the scattered field varies little (the difference between the maximum and minimum is about 0.25) from different measurement locations; that is, the scattered field is almost a constant.So the noise we add is small.If the additive noise is big corresponding to the variation, the reconstruction effect will be poor (see Figure 9).That is also why we choose two point sources illuminating together.The variation of scattered field is bigger than only one point source's situation.
In the last experiment of sound-soft boundary condition part, we explore a different prior measure.Assume the noise  ∼ N(0,  2 ),  = 0.01, and the prior density where () can be chosen as different functions, such as  1 norm function,  0 norm function, or TV norm function corresponding to different prior measures, in fact different regularization strategies.Here, we choose the  1 norm function () = ‖‖  1 and  = 7.This prior called impulse prior is corresponding to the  1 regularization.We also generate 100000 samples for simulation.The reconstruction of the cavity (Figure 10), the radial function, and the scattered field (Figure 11) are as follows.This experiment shows that Gaussian prior measure is not the only choice.The other priors may be also effective in numerical simulations.We adopt Gaussian prior measure because we know little about radial function  in advance, and the theoretical system of Gaussian prior measure is easy to build under the framework of [14].The theoretical results of other prior measures are also under consideration.
Example 2 (the sound-hard boundary condition).We truncate 6 terms of the series in (13) and assume the corresponding 13 coefficients obey Gaussian distribution N(0, 1).We take the measurement error  ∼ N(0,  2 ),  = 0.01, and generate 100000 samples of the posterior distribution.For the constraints of space, we only present the figures of reconstructed cavity, radial function, and scattered field (Figures 12  and 13).From Example 2 we can see that the Bayesian method is also effective for inverse interior scattering problem with sound-hard boundary condition.Example 2. The reconstructions of cavity, radial function, and scattered field are in Figures 14 and 15.This example shows that our method is effective in the situation of impedance boundary condition.

Conclusions
In this paper, we study the inverse interior scattering problem.Bayesian method is used to reconstruct the shape of the cavity from interior measurement.We prove the well-posedness and present some numerical examples to illustrate that our method is effective.In the future, we will consider Bayesian method for more complicated inverse scattering problem, for example, the open cavity embedded in the infinite ground plane.

Figure 1 :
Figure 1: Cavity  and the measurement location .

)
[25]f.From the above analysis of derivative   , it is easy to get that            ∞ ⩽              (B (  +   ))By mean value theorem[25]and the regularity of   and   ,       ( 1 ) −   ( 2 )