# Bayesian Method for Shape Reconstruction in the Inverse Interior Scattering Problem

Authors: Yujie Wang; Fuming Ma; Enxi Zheng
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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 D⊂R2 with Lipschitz boundary ∂D. Our aim is to reconstruct the boundary ∂D from measurements (scattered waves) taken on a closed curve C inside the cavity D (see Figure 1). The point source(s) (incident waves) are also located on curve C.Figure 1
CavityD and the measurement location C.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–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 well-posedness of Bayesian method for this problem and present some numerical results.The outline of this paper is as follows. In Section2 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.


## 2. 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 cylinderD⊂R2 be a simply connected domain with C2,α boundary ∂D. We assume the point sources and the observational points all locate on curve C inside the domain D. Then the above scattering problem reduces to finding the scattered field us which satisfies(1)Δus+k2us=0inD,(2)Bus=-Bui,on∂D,where k>0 is the wave number and ui is the incident wave given by(3)uix=Φx,d=i4H01kx-d.Here, i=-1, d is a fixed point on C, and Φ(x,d) is the fundamental solution of two dimensional Helmholtz equations. B is a linear operator defined by(4)B=I,for  sound-soft boundary condition,B=∂∂ν,for  sound-hard  boundary  condition,B=∂∂ν+iλI,for  impedance  boundary  condition.Here, λ>0 and ν is the unit outward normal to D.The direct problem is solving the above equation to get scattered fieldus for given D. The inverse problem is recovering the shape D for given measurements us. To ensure the uniqueness of interior scattering problem, we assume that k2 is not an eigenvalue for -Δ corresponding to boundary condition (2) of D. For simplification, we assume D is a starlike domain with respect to the origin. Thus, we can represent the boundary ∂D by(5)∂D≔r=rθcos⁡θ,sin⁡θ=exp⁡qθcos⁡θ,sin⁡θ,θ∈0,2π,where 0<r<R0. We set q=lnr just for the convenience of the proof in the next section. Then our model can be written as(6)y=Gq+δ,where G∶=(us(x1),…,us(xn)) (n is the number of measurements, i.e., y∈Rn) is a finite dimensional observational operator corresponding to (1) and vector y is the observational data with noise δ. We describe our prior information about q∈X, in terms of a probability measure μ0. Here, X is some function space to be chosen (see Section 3.1). And we use Bayes’ formula to calculate the posterior probability measure μy for q given y. Let π0 and πy denote the probability density functions of measures μ0 and μy. By the knowledge in statistics the noise is often a mean zero random variable and we know little about q as prior information in the interior inverse scattering problem. So we assume the noise δ~N(0,Γ) with covariance matrix Γ and the prior μ0~N(m0,C0) with mean function m0 and bounded covariance operator C0. Of course, the prior measure μ0 may be various from the different prior information of the cavity geometry q and we will present a numerical example in Section 4. According to Bayes’ formula, we can get(7)πyq∝exp⁡-12Γ-1/2y-Gq2+C0-1/2q-m02=exp⁡-12y-GqΓ2+q-m0C02,where ·Γ2∶=(·,Γ-1·) and a∝b means a is proportional to b. Then the Radon-Nikodym derivative related to the prior measure and posterior measure is(8)dμydμ0q∝exp⁡-12y-GqΓ2.In this paper, we use MCMC (Markov Chain Monte Carlo) method to describe the posterior distribution according to πy(q) and get the mean of q. For the details of this method we refer to [12, 15].


## 3. 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.


### 3.1. Well-Posedness Theorem

According to Section  4 (common structure) of [14], we should choose a space X and a prior measure μ0 such that μ0(X)=1 and Assumption 2.7 in [14] holds on X. Then we can get the well-posedness results of Bayesian method. So first we give the specific definition of prior measure on q. Assume A=-d2/dθ2 with the definition domain: (9)DA=u∈H20,2π;∫02πuθdθ=0.Because in the representation of starlike domain r(θ) 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 H2[0,2π] is the usual Sobolev space with periodic condition u(0)=u(2π); that is, (10)H20,2π≔∑n∈Zcneinθ:∑n∈Z1+n22cn2<∞.Now we assume a prior Gaussian measure on q such that q′′(θ)~N(0,A-1). Through simple calculation, we can get that the eigenvalues of A-1 are 1/n2 and the corresponding eigenfunctions are ϕ2n=cos(nθ)/π and ϕ2n-1=sin(nθ)/π. According to Karhunen-Loève expansion, we have (11)q′′θ=∑n=1∞anncos⁡nθπ+bnnsin⁡nθπ,where an and bn are i.i.d. (independent and identically distributed) sequences with a1~N(0,1) and b1~N(0,1). Clearly, this is in fact the Fourier expansion. Integrating q′′(θ) we can get q(θ). q(θ) is also a Gaussian measure because the integral operator is linear and continuous [16]. Using Lemma 6.25 in [14], we know that q′′(θ) is almost surely in C0,α for any α<1/2. In order to get q(θ), we integrate q′′(θ). The function q(θ) is not unique when its second derivative q′′(θ) is given. Here, we just use periodic form of q(θ) in (13) the same as in [10]. More specifically, we define (12)qθ=q02π-∑j=1∞ajj3π+bjj2πθ+∫0θ∫0tq′′sdsdt,where q0~N(μ,σ2) is a Gaussian random variable. Then we can integrate the Fourier expansion of q′′(θ) term by term to obtain(13)qθ=q02π-∑n=1∞ann3cos⁡nθπ+bnn3sin⁡nθπ.Here, we define the Gaussian measure on q based on its second derivative q′′. There is also an alternative way to define the Gaussian measure on q; see [10].Now we defineX=C2,α0,2π with norm ·X given in the following:(14)ψX=ψ∞+ψ′∞+ψ′′∞+supx,y∈0,2πx≠yψ′′x-ψ′′yx-yα,where ·∞ is defined as usual; that is, ψ∞=supx∈[0,2π]|ψ(x)|.Then we obtain the following theorem.Theorem 1.
Assumeq~μ0 defined by (13); then μ0(X)=1 for 0<α<1/2.Proof.
From the analysis above, we knowq′′(θ) is almost surely in C0,α for any α<1/2. So q(θ) is almost surely in C2,α[0,2π]; that is, μ0(X)=1.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 M=M(ε)∈R such that, for all q∈X,(15)GqΓ⩽exp⁡εqX2+M.
(ii) For everyt>0, there is a K=K(t)>0 such that, for all q1,q2∈X with max{q1X,q2X}<t,(16)Gq1-Gq2Γ⩽Kq1-q2X.We will prove thatG(q) satisfies Assumption 2 in the next subsection. Now we give the main theorem in our paper.Theorem 3 (well-posedness).
Let the observational operatorG satisfy Assumption 2 and q~μ0 with μ0(X)=1. Then one has the following.(i)
The posterior measureμy is absolutely continuous with respect to μ0 and has Radon-Nikodym derivative given by(17)dμydμ0q∝exp⁡-12y-GqΓ2.(ii)
The posterior measureμy is a well-defined probability measure on L2[0,2π].(iii)
The posterior measureμy is Lipschitz in the data y, with respect to the Hellinger distance; that is, there exists a constant C(t)>0 such that(18)dHellμy,μy′⩽Cty-y′Γ, for all y,y′ with max{yΓ,y′Γ<t}. Here, the Hellinger distance is defined as(19)dHellμ1,μ2≔12∫dμ1dμ0-dμ2dμ02dμ0,provided μ1 and μ2 are both absolutely continuous with μ0.The proof of Theorem3 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 μy is also continuous in y. So if we use the mean of posterior distribution to approximate the solution q(θ) of inverse interior scattering problem, it is continuous with the observational data y.


### 3.2. Properties of Observational OperatorG

In this subsection, our task is to prove Assumption2 when G(s) 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 S and K, given by(20)Sφx=2∫∂DΦx,yφydsy,x∈∂D,Kφx=2∫∂D∂Φx,y∂νyφydsy,x∈∂D,and the normal derivative operator K′, given by(21)K′φx=2∫∂D∂Φx,y∂νxφydsy,x∈∂D,where Φ(x,y) is the fundamental solution of Helmholtz equation and ν(x) is the unit outward normal.Then we look for a solutionus in the form of(22)usx=K~φx=∫∂D∂Φx,y∂νyφydsy,x∈R2∖∂D,for the Dirichlet boundary condition, or(23)usx=S~φx=∫∂DΦx,yφydsy,x∈R2∖∂D,for the Neumann boundary condition or impedance boundary condition. From the jump relations, we see that us is the solution of scattering problem (1)-(2), provided the density φ is a solution of the following integral equation:(24)φx-Kφx=2uix,x∈∂D,for the Dirichlet boundary condition, or(25)φx+K′φx=-2∂uix∂νx,x∈∂D,for the Neumann boundary condition, or(26)φx+K′φx+iλSφx=-2∂uix∂νx-2iλuix,x∈∂D,for the impedance boundary condition.From [19] we known that the operators S and K are compact operators on C2,α(∂D), and K′ is a compact operator on C1,α(∂D). If k2 is not an eigenvalue of -Δ for all possible domain D when q changes, we can get the injectivity of I-K, I+K′, and I+K′+iλS. For details we refer to [17]. In practice, the cavity is always bounded, so r(θ) is also bounded. We can choose positive constants rmin small enough and rmax big enough such that rmin<r<rmax. Then by Riesz-Fredholm theory [17], from the injectivity of the operators I-K, I+K′, and I+K′+iλS, we know that they are bijective, their inverse operators are bounded on domain D, and the bounds of the operators are only related to the constants rmin and rmax. From the analytic property of ui on ∂D, we obtain φ(x)∈C2,α(∂D) for Dirichlet boundary condition and φ(x)∈C1,α(∂D) for Neumann or impedance boundary condition. S~ maps C1,α(∂D) continuously into C2,α(D¯). As a result, us(x)∈C2,α(D¯) and the trace us(x)∂D∈C2,α(∂D) for all the three boundary conditions [19].Remark 4.
Here, for the impedance boundary condition, whenk>0 is a real number, k2 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 ak small enough to avoid all Dirichlet eigenvalues. More specifically, for r<rmax, the ball Brmax contains all possible domain D. From the Rayleigh-Faber-Krahn inequality the smallest Dirichlet eigenvalue λ1(>0) of -Δ for domain D is no less than the smallest eigenvalue λ1′(>0) of Brmax. Then when k2<λ1′, it is not an eigenvalue of -Δ for all possible domain D. 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 thatk2 is not an eigenvalue of -Δ for all possible domain D. In fact, if the cavity is filled with some special material such that Imk>0, then the interior Dirichlet and Neumann problems have at most one solution [17]. We will consider this situation in the future investigation.Now we show thatG(q) satisfies the property (i) in Assumption 2.Theorem 5.
For everyε>0, there is an M=M(ε)∈R such that, for all q∈X, (27)GqΓ⩽exp⁡εqX2+M,where |·|Γ2:=(·,Γ-1·).Proof.
From the definition ofG, we just need to prove (28)usxj⩽exp⁡εqX2+M,j=1,2,…,n.By (22) and (23), it is sufficient to show(29)∫∂DΦxj,yφydsy⩽exp⁡εqX2+M,or(30)∫∂D∂Φxj,y∂νyφydsy⩽exp⁡εqX2+M.Now we give the proof of (29); then the proof of (30) is similar. It is known that(31)∫∂DΦxj,yφydsy⩽2πΦxj,y∞φ∞1+q′2∞exp⁡q∞.Here taking Neumann boundary condition as an example, because (32)φx=-2I+K′-1∂uix∂νx,the upper bounds of φ∞ depend only on the bounds of (I+K′)-1 and ∂ui(x)/∂ν(x). From the analysis above Theorem 5, (I+K′)-1 and ∂ui(x)/∂ν(x) are bounded when the domain D is bounded. So we have φ∞≤C.
BecauseΦ(xj,y) is analytic on y∈∂D,(33)∫∂DΦxj,yφydsy⩽C1+q′2∞exp⁡q∞.
Now our goal is to bound the two factors in right hand of the above inequality.
Whenq′∞⩽1, 1+q′2∞⩽C.
Whenq′∞>1, according to Young’s inequality, (34)1+q′2∞⩽2q′∞⩽2exp⁡qX⩽2exp⁡12ε+εqX22.On the other hand, (35)exp⁡q∞⩽exp⁡qX⩽exp⁡12ε+εqX22.Then we obtain the final result: (36)GqΓ⩽exp⁡εqX2+M.The next step is to prove(37)Gq1-Gq2Γ⩽Ktq1-q2X,withmax⁡q1X,q2X<t.That is, G(q) 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–24], but the results are similar. Assume D is a bounded open domain in R2 with C2 smooth boundary ∂D. u(D) is a function defined on D satisfying (38)AuD=finD,BuD=gon∂D,where A and B are some partial differential operators. Given θ a regular (in our application regular means C2 smooth) vector field defined in R2, one denotes (39)I+θD=x+θx∣x∈D.For θ small enough, (I+θ)(D) is an open set close to D with regular boundary. (I+θ)(D) can be viewed as a variation of D. Then we define u(θ)=u((I+θ)(D)) and we call the local derivative (∂u(θ)/∂θ)(0) the domain derivative with D. Let u′=(∂uθ/∂θ)(0)τ be the local derivative of u in a direction τ; that is,(40)u′=limh→01huhτ-u0,where τ is a regular vector field.From [21] we know that, with some suitable smoothness hypotheses (e.g., C2 smooth), the derivative u′ satisfies the following equation and boundary condition:(41)∂A∂uu′=0inD,∂B∂uu′=-τ·ν∂∂νBu0-gon∂D,where ν is the normal derivative of ∂D, ∂A/∂u=(∂A(v)/∂v)(u0), and u0=u(0).In our problem, the smoothness ofD and ∂D and the regularity of u satisfy the hypotheses (3.1) to (3.9) and (3.12) in [21]. The operator A=Δ+k2I, B=B, f=0, and g=-Bui. Clearly, A and B are linear bounded operators, so the derivatives are themselves; that is, (42)∂A∂u=A,∂B∂u=B.Then we obtain the derivative u′ in our problem satisfies(43)Δu′+k2u′=0inD,(44)Bu′=-τ·∇Bus+uion∂D.Here, because of the boundary condition (2), the tangential part of ∇(B(us+ui)) is 0; that is, (45)∇Bus+ui=ν∂∂νBus-g.In our problem,D is a starlike domain. The variation of D is in fact the variation of q. So we can denote τ=exp(q)(cosθ,sinθ)q~, where q~ is the direction of variation of q. Then the boundary condition (44) changes to (46)Bu′=-∂Bus+ui∂erexp⁡qq~,where er=(cosθ,sinθ).Now we show thatG(q) satisfies property (ii) in Assumption 2.Theorem 6.
For everyt>0, there is a K=K(t)>0 such that, for all q1,q2∈X with max{q1X,q2X}<t, (47)Gq1-Gq2Γ⩽Kq1-q2X.Proof.
From the above analysis of derivativeu′, it is easy to get that (48)u′∞⩽C∂Bus+ui∂er∞exp⁡q∞q~∞.By mean value theorem [25] and the regularity of us and ui, (49)usq1-usq2∞⩽u′∞q1-q2∞⩽Cexp⁡max⁡q1∞,q2∞q1-q2∞.Then from the definition of G, we obtain (50)Gq1-Gq2Γ⩽Cexp⁡max⁡q1∞,q2∞q1-q2∞⩽Cexp⁡max⁡q1X,q2Xq1-q2X.


## 3.1. Well-Posedness Theorem

According to Section  4 (common structure) of [14], we should choose a space X and a prior measure μ0 such that μ0(X)=1 and Assumption 2.7 in [14] holds on X. Then we can get the well-posedness results of Bayesian method. So first we give the specific definition of prior measure on q. Assume A=-d2/dθ2 with the definition domain: (9)DA=u∈H20,2π;∫02πuθdθ=0.Because in the representation of starlike domain r(θ) 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 H2[0,2π] is the usual Sobolev space with periodic condition u(0)=u(2π); that is, (10)H20,2π≔∑n∈Zcneinθ:∑n∈Z1+n22cn2<∞.Now we assume a prior Gaussian measure on q such that q′′(θ)~N(0,A-1). Through simple calculation, we can get that the eigenvalues of A-1 are 1/n2 and the corresponding eigenfunctions are ϕ2n=cos(nθ)/π and ϕ2n-1=sin(nθ)/π. According to Karhunen-Loève expansion, we have (11)q′′θ=∑n=1∞anncos⁡nθπ+bnnsin⁡nθπ,where an and bn are i.i.d. (independent and identically distributed) sequences with a1~N(0,1) and b1~N(0,1). Clearly, this is in fact the Fourier expansion. Integrating q′′(θ) we can get q(θ). q(θ) is also a Gaussian measure because the integral operator is linear and continuous [16]. Using Lemma 6.25 in [14], we know that q′′(θ) is almost surely in C0,α for any α<1/2. In order to get q(θ), we integrate q′′(θ). The function q(θ) is not unique when its second derivative q′′(θ) is given. Here, we just use periodic form of q(θ) in (13) the same as in [10]. More specifically, we define (12)qθ=q02π-∑j=1∞ajj3π+bjj2πθ+∫0θ∫0tq′′sdsdt,where q0~N(μ,σ2) is a Gaussian random variable. Then we can integrate the Fourier expansion of q′′(θ) term by term to obtain(13)qθ=q02π-∑n=1∞ann3cos⁡nθπ+bnn3sin⁡nθπ.Here, we define the Gaussian measure on q based on its second derivative q′′. There is also an alternative way to define the Gaussian measure on q; see [10].Now we defineX=C2,α0,2π with norm ·X given in the following:(14)ψX=ψ∞+ψ′∞+ψ′′∞+supx,y∈0,2πx≠yψ′′x-ψ′′yx-yα,where ·∞ is defined as usual; that is, ψ∞=supx∈[0,2π]|ψ(x)|.Then we obtain the following theorem.Theorem 1.
Assumeq~μ0 defined by (13); then μ0(X)=1 for 0<α<1/2.Proof.
From the analysis above, we knowq′′(θ) is almost surely in C0,α for any α<1/2. So q(θ) is almost surely in C2,α[0,2π]; that is, μ0(X)=1.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 M=M(ε)∈R such that, for all q∈X,(15)GqΓ⩽exp⁡εqX2+M.
(ii) For everyt>0, there is a K=K(t)>0 such that, for all q1,q2∈X with max{q1X,q2X}<t,(16)Gq1-Gq2Γ⩽Kq1-q2X.We will prove thatG(q) satisfies Assumption 2 in the next subsection. Now we give the main theorem in our paper.Theorem 3 (well-posedness).
Let the observational operatorG satisfy Assumption 2 and q~μ0 with μ0(X)=1. Then one has the following.(i)
The posterior measureμy is absolutely continuous with respect to μ0 and has Radon-Nikodym derivative given by(17)dμydμ0q∝exp⁡-12y-GqΓ2.(ii)
The posterior measureμy is a well-defined probability measure on L2[0,2π].(iii)
The posterior measureμy is Lipschitz in the data y, with respect to the Hellinger distance; that is, there exists a constant C(t)>0 such that(18)dHellμy,μy′⩽Cty-y′Γ, for all y,y′ with max{yΓ,y′Γ<t}. Here, the Hellinger distance is defined as(19)dHellμ1,μ2≔12∫dμ1dμ0-dμ2dμ02dμ0,provided μ1 and μ2 are both absolutely continuous with μ0.The proof of Theorem3 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 μy is also continuous in y. So if we use the mean of posterior distribution to approximate the solution q(θ) of inverse interior scattering problem, it is continuous with the observational data y.


## 3.2. Properties of Observational OperatorG

In this subsection, our task is to prove Assumption2 when G(s) 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 S and K, given by(20)Sφx=2∫∂DΦx,yφydsy,x∈∂D,Kφx=2∫∂D∂Φx,y∂νyφydsy,x∈∂D,and the normal derivative operator K′, given by(21)K′φx=2∫∂D∂Φx,y∂νxφydsy,x∈∂D,where Φ(x,y) is the fundamental solution of Helmholtz equation and ν(x) is the unit outward normal.Then we look for a solutionus in the form of(22)usx=K~φx=∫∂D∂Φx,y∂νyφydsy,x∈R2∖∂D,for the Dirichlet boundary condition, or(23)usx=S~φx=∫∂DΦx,yφydsy,x∈R2∖∂D,for the Neumann boundary condition or impedance boundary condition. From the jump relations, we see that us is the solution of scattering problem (1)-(2), provided the density φ is a solution of the following integral equation:(24)φx-Kφx=2uix,x∈∂D,for the Dirichlet boundary condition, or(25)φx+K′φx=-2∂uix∂νx,x∈∂D,for the Neumann boundary condition, or(26)φx+K′φx+iλSφx=-2∂uix∂νx-2iλuix,x∈∂D,for the impedance boundary condition.From [19] we known that the operators S and K are compact operators on C2,α(∂D), and K′ is a compact operator on C1,α(∂D). If k2 is not an eigenvalue of -Δ for all possible domain D when q changes, we can get the injectivity of I-K, I+K′, and I+K′+iλS. For details we refer to [17]. In practice, the cavity is always bounded, so r(θ) is also bounded. We can choose positive constants rmin small enough and rmax big enough such that rmin<r<rmax. Then by Riesz-Fredholm theory [17], from the injectivity of the operators I-K, I+K′, and I+K′+iλS, we know that they are bijective, their inverse operators are bounded on domain D, and the bounds of the operators are only related to the constants rmin and rmax. From the analytic property of ui on ∂D, we obtain φ(x)∈C2,α(∂D) for Dirichlet boundary condition and φ(x)∈C1,α(∂D) for Neumann or impedance boundary condition. S~ maps C1,α(∂D) continuously into C2,α(D¯). As a result, us(x)∈C2,α(D¯) and the trace us(x)∂D∈C2,α(∂D) for all the three boundary conditions [19].Remark 4.
Here, for the impedance boundary condition, whenk>0 is a real number, k2 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 ak small enough to avoid all Dirichlet eigenvalues. More specifically, for r<rmax, the ball Brmax contains all possible domain D. From the Rayleigh-Faber-Krahn inequality the smallest Dirichlet eigenvalue λ1(>0) of -Δ for domain D is no less than the smallest eigenvalue λ1′(>0) of Brmax. Then when k2<λ1′, it is not an eigenvalue of -Δ for all possible domain D. 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 thatk2 is not an eigenvalue of -Δ for all possible domain D. In fact, if the cavity is filled with some special material such that Imk>0, then the interior Dirichlet and Neumann problems have at most one solution [17]. We will consider this situation in the future investigation.Now we show thatG(q) satisfies the property (i) in Assumption 2.Theorem 5.
For everyε>0, there is an M=M(ε)∈R such that, for all q∈X, (27)GqΓ⩽exp⁡εqX2+M,where |·|Γ2:=(·,Γ-1·).Proof.
From the definition ofG, we just need to prove (28)usxj⩽exp⁡εqX2+M,j=1,2,…,n.By (22) and (23), it is sufficient to show(29)∫∂DΦxj,yφydsy⩽exp⁡εqX2+M,or(30)∫∂D∂Φxj,y∂νyφydsy⩽exp⁡εqX2+M.Now we give the proof of (29); then the proof of (30) is similar. It is known that(31)∫∂DΦxj,yφydsy⩽2πΦxj,y∞φ∞1+q′2∞exp⁡q∞.Here taking Neumann boundary condition as an example, because (32)φx=-2I+K′-1∂uix∂νx,the upper bounds of φ∞ depend only on the bounds of (I+K′)-1 and ∂ui(x)/∂ν(x). From the analysis above Theorem 5, (I+K′)-1 and ∂ui(x)/∂ν(x) are bounded when the domain D is bounded. So we have φ∞≤C.
BecauseΦ(xj,y) is analytic on y∈∂D,(33)∫∂DΦxj,yφydsy⩽C1+q′2∞exp⁡q∞.
Now our goal is to bound the two factors in right hand of the above inequality.
Whenq′∞⩽1, 1+q′2∞⩽C.
Whenq′∞>1, according to Young’s inequality, (34)1+q′2∞⩽2q′∞⩽2exp⁡qX⩽2exp⁡12ε+εqX22.On the other hand, (35)exp⁡q∞⩽exp⁡qX⩽exp⁡12ε+εqX22.Then we obtain the final result: (36)GqΓ⩽exp⁡εqX2+M.The next step is to prove(37)Gq1-Gq2Γ⩽Ktq1-q2X,withmax⁡q1X,q2X<t.That is, G(q) 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–24], but the results are similar. Assume D is a bounded open domain in R2 with C2 smooth boundary ∂D. u(D) is a function defined on D satisfying (38)AuD=finD,BuD=gon∂D,where A and B are some partial differential operators. Given θ a regular (in our application regular means C2 smooth) vector field defined in R2, one denotes (39)I+θD=x+θx∣x∈D.For θ small enough, (I+θ)(D) is an open set close to D with regular boundary. (I+θ)(D) can be viewed as a variation of D. Then we define u(θ)=u((I+θ)(D)) and we call the local derivative (∂u(θ)/∂θ)(0) the domain derivative with D. Let u′=(∂uθ/∂θ)(0)τ be the local derivative of u in a direction τ; that is,(40)u′=limh→01huhτ-u0,where τ is a regular vector field.From [21] we know that, with some suitable smoothness hypotheses (e.g., C2 smooth), the derivative u′ satisfies the following equation and boundary condition:(41)∂A∂uu′=0inD,∂B∂uu′=-τ·ν∂∂νBu0-gon∂D,where ν is the normal derivative of ∂D, ∂A/∂u=(∂A(v)/∂v)(u0), and u0=u(0).In our problem, the smoothness ofD and ∂D and the regularity of u satisfy the hypotheses (3.1) to (3.9) and (3.12) in [21]. The operator A=Δ+k2I, B=B, f=0, and g=-Bui. Clearly, A and B are linear bounded operators, so the derivatives are themselves; that is, (42)∂A∂u=A,∂B∂u=B.Then we obtain the derivative u′ in our problem satisfies(43)Δu′+k2u′=0inD,(44)Bu′=-τ·∇Bus+uion∂D.Here, because of the boundary condition (2), the tangential part of ∇(B(us+ui)) is 0; that is, (45)∇Bus+ui=ν∂∂νBus-g.In our problem,D is a starlike domain. The variation of D is in fact the variation of q. So we can denote τ=exp(q)(cosθ,sinθ)q~, where q~ is the direction of variation of q. Then the boundary condition (44) changes to (46)Bu′=-∂Bus+ui∂erexp⁡qq~,where er=(cosθ,sinθ).Now we show thatG(q) satisfies property (ii) in Assumption 2.Theorem 6.
For everyt>0, there is a K=K(t)>0 such that, for all q1,q2∈X with max{q1X,q2X}<t, (47)Gq1-Gq2Γ⩽Kq1-q2X.Proof.
From the above analysis of derivativeu′, it is easy to get that (48)u′∞⩽C∂Bus+ui∂er∞exp⁡q∞q~∞.By mean value theorem [25] and the regularity of us and ui, (49)usq1-usq2∞⩽u′∞q1-q2∞⩽Cexp⁡max⁡q1∞,q2∞q1-q2∞.Then from the definition of G, we obtain (50)Gq1-Gq2Γ⩽Cexp⁡max⁡q1∞,q2∞q1-q2∞⩽Cexp⁡max⁡q1X,q2Xq1-q2X.


## 4. 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 byr(t)=0.44cos2t+sin2t. The interior curve C is a circle defined by C={ρ(t)∣ρ(t)=0.3(cost,sint),t∈[0,2π]}. We take the source points d1=(0,0.3) and d2=(0,-0.3) on C. The wave number k=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(q)) on curve C. 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 q0 and a1 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 q0 and a1. 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 r(t) and the scattered field us are displayed in Figure 5, where the x-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 Figures6 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 (51)π0q∝exp-βHq,where H(q) can be chosen as different functions, such as L1 norm function, L0 norm function, or TV norm function corresponding to different prior measures, in fact different regularization strategies. Here, we choose the L1 norm function H(q)=qL1 and β=7. This prior called impulse prior is corresponding to the L1 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 r 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.Figure 2
Posterior distributions ofq0 and a1.
(a)
(b)Figure 3
Markov chains ofq0 and a1.
(a)
(b)Figure 4
Reconstruction of the cavity (σ=0.005).Figure 5
Reconstruction of the radial function (a) and scattered field (b) (σ=0.005).
(a)
(b)Figure 6
Reconstruction of the cavity (σ=0.01).Figure 7
Reconstruction of the radial function (a) and scattered field (b) (σ=0.01).
(a)
(b)Figure 8
Reconstruction of the cavity (σ=0.05).Figure 9
Reconstruction of the radial function (a) and scattered field (b) (σ=0.05).
(a)
(b)Figure 10
Reconstruction of the cavity (L1 prior, σ=0.01).Figure 11
Reconstruction of the radial function (a) and scattered field (b) (L1 prior, σ=0.01).
(a)
(b)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.Figure 12
Reconstruction of the cavity (sound-hard boundary condition,σ=0.01).Figure 13
Reconstruction of the radial function (a) and scattered field (b) (sound-hard boundary condition,σ=0.01).
(a)
(b)Example 3 (the impedance boundary condition).
We chooseλ=1 in this example and other parameters are the same as in 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.Figure 14
Reconstruction of the cavity (impedance boundary condition,σ=0.01).Figure 15
Reconstruction of the radial function (a) and scattered field (b) (impedance boundary condition,σ=0.01).
(a)
(b)


## 5. 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.



---

*Source: 935294-2015-05-28.xml*