# Two-Step Extended Sampling Method for the Inverse Acoustic Source Problem

Authors: Jiyu Sun; Yuhui Han
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Recently, a new method, called the extended sampling method (ESM), was proposed for the inverse scattering problems. Similar to the classical linear sampling method (LSM), the ESM is simple to implement and fast. Compared to the LSM which uses full-aperture scattering data, the ESM only uses the scattering data of one incident wave. In this paper, we generalize the ESM for the inverse acoustic source problems. We show that the indicator function of ESM, which is defined using the approximated solutions of some linear ill-posed integral equations, is small when the support of the source is contained in the sampling disc and is large when the source is outside. This behavior is similar to the ESM for the inverse scattering problem. Numerical examples are presented to show the effectiveness of the method.


---

## Body

## 1. Introduction

Recently, a new method, called the extended sampling method (ESM), was proposed for the inverse scattering problems using a single incident wave at a fixed frequency [1]. ESM is based on an idea similar to that of the linear sampling method, which has been studied extensively in the literature [2, 3]. Although the ESM can only reconstruct the location and approximate the size of the target, it uses much less data than the classical linear sampling method, which needs full-aperture data (far-field pattern of all incident directions and all observation directions). This property makes the ESM suitable for the inverse source problems.Inverse source problems (ISP) have many applications such as pollution source detection, optical tomography, and sound source localization [4–8]. It is known that the inverse source problem using a single frequency far-field data does not have a unique solution [6, 9]. In this paper, we consider the inverse problem to determine the location and size of an acoustic source at a fixed frequency. The ESM is generalized for the inverse acoustic source problems. We show that the indicator function of the proposed method, which is defined using the approximated solutions of some linear ill-posed integral equations, is small when the support of the source is contained in the sampling disc and is large when the source is outside. The behavior of the indicator function is similar to the ESM for the inverse scattering problem. Numerical examples are presented to validate the effectiveness of the method.The ESM was first proposed for the inverse scattering problem to reconstruct the obstacles [1]. To treat the inverse source problem, we introduced an artificial intermediate scattering problem, whose far-field pattern is the same as the original source problem. This enables us to generalize the ESM for the inverse source problems with little change in the numerical algorithm.The rest of this paper is arranged as follows. In Section2, we introduce the acoustic source problems and give a brief presentation of the ESM for the inverse scattering problems. Section 3 contains the generalization of the ESM for the inverse source problem and its mathematical justification. In particular, we show that the indicator function of the ESM is small when the support of the source is contained in the sampling disc and is large when the source is outside. Section 4 contains several validating numerical examples in two dimensions. We make some conclusions in Section 5.


## 2. Preliminaries

In this section, we describe the acoustic source problem and the associated inverse problem of interests. Then, we briefly introduce the ESM for the inverse acoustic scattering problem.


### 2.1. Acoustic Source Problem

LetD⊂ℝ2 be a bounded domain and k>0 be the wave number. Let fx∈L2D represent the source with compact support D. The acoustic source problem is to find a function u∈Hloc1ℝ2 satisfying the Helmholtz equation(1)Δu+k2u=fx,inℝ2,and the Sommerfeld radiation condition(2)limr⟶∞r1/2∂u∂r−iku=0,r=x.The above problem has a solutionu given by(3)ux=∫DΦx,yfydy,x∈ℝ2,where(4)Φx,y=i4H01kx−y,is the fundamental solution of the Helmholtz equation. Here, H01 is the Hankel function of the first kind and order 0 given by(5)H01=J0+iY0,where the real part J0 is the Bessel function of order 0 and the imaginary part Y0 is the Neumann function of order 0 (see page 73 of [3]).Since the solutionu to (1) is radiating, it has the following asymptotic behavior (see, e.g., Theorem 2.6 in [3]):(6)ux=eiπ/48πkeikxx1/2u∞x^+Ox−3/2,r=x⟶∞,where x^=x/x∈S≔x^∈ℝ2x^=1 and u∞x^ is called the far-field pattern of u. From the asymptotic behavior of the Hankel function H01 and (3), we have that(7)u∞x^=∫De−ikx^⋅yfydy,x^∈S.Theinverse source problem (ISP) considered in this paper is to find the location and size of D from the far-field pattern u∞x^.


### 2.2. Extended Sampling Method

The ESM was first proposed in [1] to obtain the location and size of an obstacle using the scattering data of one incident plane wave. We first give a quick introduction for the scattering problem of a sound-soft obstacle. Note that the method works for other types of obstacles.We still useD⊂ℝ2 to denote a bounded domain occupied by the sound-soft obstacle. Let uix≔eikx⋅d,x,d∈ℝ2 be the incident plane wave, where d=1 is the direction. The direct exterior scattering problem is to find the scattered field us such that(8)Δus+k2us=0inℝ2\D¯,us=−uion∂D,limr⟶∞r∂us∂r−ikus=0.Since the scattered fieldus is radiating, it has an asymptotic expansion(9)usx=eiπ/48kπeikrru∞x^+O1r,asr≔x⟶∞,uniformly in all directions x^=x/x. The function u∞x^ is the far-field pattern of us due to the incident field ui.When the obstacle is a disc centered at the origin with radiusR, denoted by B, the far-field pattern can be written down exactly (see, e.g., Chp. 3 of [3]):(10)u∞Bx^;d=−e−iπ/42πkJ0kRH01kR+2∑n=1∞JnkRHn1kRcosnθ,x^∈S,where Jn is the Bessel function, Hn1 is the Hankel function of the first kind of order n, and θ=∠x^,d is the angle between x^ and d. If the disc is shifted to z, i.e.,(11)Bz≔x+z;x∈B,z∈ℝ2,the far-field pattern is simply(12)u∞Bzx^;d=eikz⋅d−x^u∞Bx^,d,x^∈S.The inverse scattering problem is to find the location and size ofD from the far-field pattern u∞x^ due to a single incident plane wave. For a point z∈ℝ2, the ESM defines a far-field operator ℱz:L2S⟶L2S using the far-field pattern of the disc Bz due to the incident plane waves of all directions d∈S:(13)ℱzgx^=∫Su∞Bzx^,dgddsd,x^∈S,where s⋅ indicates the surface integral.Letusx and u∞x^ be the scattered field and far-field pattern of the scatterer D due to one incident wave, respectively. Using the far-field operator ℱz, we set up a far-field equation(14)ℱzgx^=u∞x^,x^∈S.Suppose thatΩ is a domain having D inside. One can generate a set of sampling points S for Ω. For z∈S, the following result is proved in [1]. Let Bz be a sound-soft disc centered at z with radius R and assume that k2 is not a Dirichlet eigenvalue for Bz. For the far-field equation (14), if D⊂Bz, for a given ε>0, there exists a function gzε∈L2S such that(15)∫Su∞Bzx^,dgzεddsd−u∞x^L2S<ε,and the Herglotz wave function(16)vgzεx≔∫Seikx⋅dgzεddsd,x∈Bz,converges to the solution w∈H1Bz of the Helmholtz equation with w=−us on ∂Bz as ε⟶0. Otherwise, if D∩Bz=∅, every gzε∈L2S that satisfies (15) for a given ε>0 is such that(17)limε⟶0vgzεH1Bz=∞.Remark 1.
The valuek2 is called a Dirichlet eigenvalue for a bounded domain B if there exists at least one nontrivial solution u to the following problem:(18)Δu+k2u=0,inB,u=0,on∂B.
Using (12), u∞Bz⋅,d can be computed easily. The far-field pattern u∞x^ is the measured data due to an incident plane wave. For each sampling point z∈S, (14) can be solved by some regularization scheme, say Tikhonov regularization. From the above discussion, the approximate solution gzαL2S should be relatively large when D is not inside Bz and relatively small when D is inside Bz. Then, an approximation of the location and support of D can be obtained by plotting gzαL2S for all sampling points z∈S.Remark 2.
Equation (14) is similar to that of the linear sampling method by Colton and Kirsch [2]. The far-field equation of the linear sampling method is(19)∫Su∞x^;dgzddsd=Φ∞x^;z,x^∈S,where u∞x^;d,x^,d∈S is the full-aperture far-field pattern of the obstacle D due to the incident plane wave ui=expikx⋅d. For a sampling point z∈Ω which contains D, using some appropriate regularization scheme, one obtains an approximate solution gzα. In general, the norm gzα is smaller for z∈D and larger for z∉D. Since u∞x^;d is the kernel in (19), the linear sampling method cannot directly process the far-field pattern due to a single incident plane wave.


## 2.1. Acoustic Source Problem

LetD⊂ℝ2 be a bounded domain and k>0 be the wave number. Let fx∈L2D represent the source with compact support D. The acoustic source problem is to find a function u∈Hloc1ℝ2 satisfying the Helmholtz equation(1)Δu+k2u=fx,inℝ2,and the Sommerfeld radiation condition(2)limr⟶∞r1/2∂u∂r−iku=0,r=x.The above problem has a solutionu given by(3)ux=∫DΦx,yfydy,x∈ℝ2,where(4)Φx,y=i4H01kx−y,is the fundamental solution of the Helmholtz equation. Here, H01 is the Hankel function of the first kind and order 0 given by(5)H01=J0+iY0,where the real part J0 is the Bessel function of order 0 and the imaginary part Y0 is the Neumann function of order 0 (see page 73 of [3]).Since the solutionu to (1) is radiating, it has the following asymptotic behavior (see, e.g., Theorem 2.6 in [3]):(6)ux=eiπ/48πkeikxx1/2u∞x^+Ox−3/2,r=x⟶∞,where x^=x/x∈S≔x^∈ℝ2x^=1 and u∞x^ is called the far-field pattern of u. From the asymptotic behavior of the Hankel function H01 and (3), we have that(7)u∞x^=∫De−ikx^⋅yfydy,x^∈S.Theinverse source problem (ISP) considered in this paper is to find the location and size of D from the far-field pattern u∞x^.


## 2.2. Extended Sampling Method

The ESM was first proposed in [1] to obtain the location and size of an obstacle using the scattering data of one incident plane wave. We first give a quick introduction for the scattering problem of a sound-soft obstacle. Note that the method works for other types of obstacles.We still useD⊂ℝ2 to denote a bounded domain occupied by the sound-soft obstacle. Let uix≔eikx⋅d,x,d∈ℝ2 be the incident plane wave, where d=1 is the direction. The direct exterior scattering problem is to find the scattered field us such that(8)Δus+k2us=0inℝ2\D¯,us=−uion∂D,limr⟶∞r∂us∂r−ikus=0.Since the scattered fieldus is radiating, it has an asymptotic expansion(9)usx=eiπ/48kπeikrru∞x^+O1r,asr≔x⟶∞,uniformly in all directions x^=x/x. The function u∞x^ is the far-field pattern of us due to the incident field ui.When the obstacle is a disc centered at the origin with radiusR, denoted by B, the far-field pattern can be written down exactly (see, e.g., Chp. 3 of [3]):(10)u∞Bx^;d=−e−iπ/42πkJ0kRH01kR+2∑n=1∞JnkRHn1kRcosnθ,x^∈S,where Jn is the Bessel function, Hn1 is the Hankel function of the first kind of order n, and θ=∠x^,d is the angle between x^ and d. If the disc is shifted to z, i.e.,(11)Bz≔x+z;x∈B,z∈ℝ2,the far-field pattern is simply(12)u∞Bzx^;d=eikz⋅d−x^u∞Bx^,d,x^∈S.The inverse scattering problem is to find the location and size ofD from the far-field pattern u∞x^ due to a single incident plane wave. For a point z∈ℝ2, the ESM defines a far-field operator ℱz:L2S⟶L2S using the far-field pattern of the disc Bz due to the incident plane waves of all directions d∈S:(13)ℱzgx^=∫Su∞Bzx^,dgddsd,x^∈S,where s⋅ indicates the surface integral.Letusx and u∞x^ be the scattered field and far-field pattern of the scatterer D due to one incident wave, respectively. Using the far-field operator ℱz, we set up a far-field equation(14)ℱzgx^=u∞x^,x^∈S.Suppose thatΩ is a domain having D inside. One can generate a set of sampling points S for Ω. For z∈S, the following result is proved in [1]. Let Bz be a sound-soft disc centered at z with radius R and assume that k2 is not a Dirichlet eigenvalue for Bz. For the far-field equation (14), if D⊂Bz, for a given ε>0, there exists a function gzε∈L2S such that(15)∫Su∞Bzx^,dgzεddsd−u∞x^L2S<ε,and the Herglotz wave function(16)vgzεx≔∫Seikx⋅dgzεddsd,x∈Bz,converges to the solution w∈H1Bz of the Helmholtz equation with w=−us on ∂Bz as ε⟶0. Otherwise, if D∩Bz=∅, every gzε∈L2S that satisfies (15) for a given ε>0 is such that(17)limε⟶0vgzεH1Bz=∞.Remark 1.
The valuek2 is called a Dirichlet eigenvalue for a bounded domain B if there exists at least one nontrivial solution u to the following problem:(18)Δu+k2u=0,inB,u=0,on∂B.
Using (12), u∞Bz⋅,d can be computed easily. The far-field pattern u∞x^ is the measured data due to an incident plane wave. For each sampling point z∈S, (14) can be solved by some regularization scheme, say Tikhonov regularization. From the above discussion, the approximate solution gzαL2S should be relatively large when D is not inside Bz and relatively small when D is inside Bz. Then, an approximation of the location and support of D can be obtained by plotting gzαL2S for all sampling points z∈S.Remark 2.
Equation (14) is similar to that of the linear sampling method by Colton and Kirsch [2]. The far-field equation of the linear sampling method is(19)∫Su∞x^;dgzddsd=Φ∞x^;z,x^∈S,where u∞x^;d,x^,d∈S is the full-aperture far-field pattern of the obstacle D due to the incident plane wave ui=expikx⋅d. For a sampling point z∈Ω which contains D, using some appropriate regularization scheme, one obtains an approximate solution gzα. In general, the norm gzα is smaller for z∈D and larger for z∉D. Since u∞x^;d is the kernel in (19), the linear sampling method cannot directly process the far-field pattern due to a single incident plane wave.


## 3. Two-Step ESM for the ISP

Now we generalize the ESM for the inverse source problems. Letu∞x^ be the far-field pattern for the inverse source problem defined in (7). Since the solution of (1) is radiating, we can use the same methodology as for the inverse scattering problems. The far-field equation for the inverse source problem keeps the same as (9), i.e.,(20)∫Su∞Bzx^,dgzddsd=u∞x^,x^∈S.The only difference is that, in (20), the far-field pattern is for the acoustic source problem defined in Section 2.1.The following theorem holds for the inverse source problem.Theorem 1.
LetBz be a sound-soft disc centered at z with radius R. Assume that the radius R is large enough such that D¯⊂Bz for some z∈ℝ2. Furthermore, assume that k2 is not a Dirichlet eigenvalue for Bz. Let u∞x^ be the far-field pattern for the acoustic source problem defined in (7). Then, the following properties hold for the far-field equation (14):Case 1.
IfD¯⊂Bz, for a given ε>0, there exists a function gzε∈L2S such that(21)∫Su∞Bzx^,dgzεddsd−u∞x^L2S<ε,and the Herglotz wave function(22)vgzεx≔∫Seikx⋅dgzεddsd,x∈Bz,converges to the solution w∈H1Bz of the Helmholtz equation with w=u on ∂Bz as ε⟶0.Case 2.
IfD¯∩B¯z=∅, every gzε∈L2S that satisfies (15) for a given ε>0 is such that(23)limε⟶0vgzεH1Bz=∞.Proof.
It is easy to see that, under the conditions of Cases1 and 2, one can find a domain D0 such that D⊂D0. Define(24)d1=maxx∈∂D,y∈∂D0x−y,d2=minx∈∂D,y∈∂D0x−y.
SinceD¯⊂Bz, there exist a positive number σ small enough such that 0<d2<d1<σ. Consequently, we can further require that D0⊂Bz for some z∈ℝ2.
Due to the well-posedness of the acoustic source problem, there exists a radiating solutionu for (1). We define an auxiliary scattering problem for a sound-soft obstacle of finding the scattered field us such that(25)Δu0s+k2u0s=0inℝ2\D¯0,u0s=u∂D0on∂D0,limr⟶∞r∂u0s∂r−iku0s=0.
Then, there exists a radiating solutionus to the above equation. Furthermore, us coincides with u in ℝ2\D¯0. By introducing D0 and the related scattering problem, one actually deals with an inverse scattering problem instead of the inverse source problem. The rest of the proof is then based on the idea of Theorem 3.1 of [1]. For completeness, we present the details as follows.
Recall that the Herglotz wave operatorℋ:L2S⟶H1/2∂Bz is defined as(26)ℋgx≔∫Seikx⋅dgddsd,x∈∂Bz.
It is well known thatℋ is injective and has a dense range provided k2 is not a Dirichlet eigenvalue for the negative Laplacian for Bz [3].
Define the near-to-far-field operator by(27)N:H1/2∂Bz⟶L2S,which maps the boundary value of a radiating solution u∈Hloc1ℝ2\B¯z to its far-field pattern u∞. Then, the operator N is bounded and injective and has a dense range [3].
Consider the case ofD¯⊂Bz. For any ε, there exists gzε such that(28)ℋgzε+usL2∂Bz≤εN.
Then,(29)N−ℋgzε=∫Su∞Bz⋅,dgzεddsd.
SinceD0⊂Bz, we also have(30)NUs=u∞.
Subtracting (29) from (30), and using (28), we have gzε such that(31)∫Su∞Bzx^,dgzεddsd−u∞x^L2S=N−ℋgzε−NusL2S≤N⋅ℋgzε+usL2S≤ε,for any ε. Taking ε⟶0 in (28), the Herglotz wave function(32)vgzεx≔∫Seikx⋅dgzεddsd,x∈Bz,converges to the unique solution w∈H1Bz of the following problem(33)Δw+k2w=0,inBz,w=−us,on∂Bz,which proves Case 1.
For the case ofD¯∩B¯z=∅, there exists gzε∈L2S such that(34)∫Su∞Bzx^,dgzεddsd−u∞x^L2S<ε,for every ε>0. Assume to the contrary that there exists a sequence vgzεn such that vgzεnH1Bz remains bounded as εn⟶0,n⟶∞. Without loss of generality, we assume that vgzεn converges to vgz∈H1Bz weakly as n⟶∞, where(35)vgzx=∫Seikx⋅dgzddsd,x∈Bz.
Letvs∈Hloc1ℝ2\B¯z be the unique solution of the following problem:(36)Δvs+k2vs=0inℝ2\B¯z,vs=−vgzon∂Bz,limr⟶∞r∂vs∂r−ikvs=0.
Its far-field pattern is(37)v∞=∫Su∞Bzx,dgzddsd.
From (34), as ε⟶0, we have(38)∫Su∞Bzx,dgzddsd=U∞.
Consequently,(39)v∞=u∞.
Then by Rellich’s lemma (see Lemma 2.12 in [3]), the scattered fields coincide in ℝ2\B¯z∪D¯0 and(40)u˜s≔vs=us,inℝ2\B¯z∪D¯0.
We have thatvs has an extension into ℝ2\Bz and us has an extension into ℝ2\D0. Since D0∩Bz=∅, U˜s can be extended from ℝ2\B¯z∪D0¯ into all of ℝ2; that is, u˜s is an entire solution to the Helmholtz equation. Since u˜s also satisfies the radiation condition, it must vanish identically in all of ℝ2. This leads to a contradiction since u˜s does not vanish identically.Remark 3.
The case ofD∩Bz≠∅ and D⊄Bz is not covered by the above theorem. In fact, the behavior of the regularized solution gz can be analyzed if one knows D. However, this is irrelevant for the inverse problem since D is always unknown.
The above theorem provides a guidance to determine the location of the source. However, the choice ofR is rather arbitrary. To find both the location and rough size of the source, we propose a two-step ESM as follows.
Two-Step ESM for the ISP(A)
Find the location ofD.(i)
Assume that a domainΩ containing D is known a priori.(ii)
Choose a large enough radiusR for Bz and generate a set S of sampling points for Ω.(iii)
Computeu∞Bzx^,d for all x^∈S and d∈S.(iv)
For each sampling pointz∈S,(a) set up a discrete version of (20);(b) use the Tikhonov regularization to compute an approximategzα to (20) (α is the regularization parameter).(v)
Take the global minimum pointz∗∈S for gzL2α as the location for D.(B)
Find the size ofD.(i)
Generate a decreasing sequence ofRj,j=1,2,…,J such that R>R1>⋯>RJ.(ii)
Forj=1,2,…,J,(a) generateSj with the distance between sampling points being roughly Rj;(b) carry out Step A usingSj and a sampling disc with radius Rj;(c) find the minimum pointzj∗∈Sj and set Dj≔BzjRj. If zj∉Dj−1, stop.(iii)
Choosezj−1 and Bzj−1Rj−1 to be the location and approximate support of D, respectively.Remark 4.
The requirement thatk2 is not a Dirichlet eigenvalue of Bz is not essential. One can avoid this by choosing a slightly different R.


## 4. Numerical Examples

In this section, we present some numerical examples to demonstrate the effectiveness of the proposed method. The synthetic data of the direct source problems are computed using the integral formula (7). Two wavenumbers k=1 and k=6 are selected for tests. We first generate a triangular mesh T for D with mesh size h≈0.01 (see Figure 1).Figure 1
Sample meshes for two test domains.
(a)
(b)The far-field pattern is approximated by(41)u∞x^≈∑T∈Te−ikx^⋅yTfyTT,where T∈T is a triangle, yT is the center of T, and T is the area of T. The interval 0,2π is uniformly divided into 40 intervals, and let θj=2jπ/40,j=0,…,39. For all examples, assume that u∞x^j,x^j=cosθj,sinθjT are the measured data. Note that the quadrature rule leads to about 3% error for the far-field pattern, which can be viewed as the noise.Let the sampling domain be given byΩ≔−8,8×−8,8. The interval −8,8 is equally divided into 160 intervals and we end up with 161×161 sampling points uniformly distributed in Ω. We denote by S the set of all sampling points. For each point z∈S, compute(42)u∞Bzx^i;dj,x^i=cosθi,sinθi,dj=cosθj,sinθjT,i,j=0,1,…,39,using (12). Simply applying the trapezoidal rule to (20), one gets the following linear system:(43)Ug=u,where(44)Uij=π20u∞Bzx^i;dj,i,j=0,1,…,39,g=gzd1,…gzd39T,j=0,1,…,39,u=u∞dj,j=0,1,…,39.Then, the Tikhonov regularization is used to solve an approximated solutiongzα of (43) with a fixed parameter α=10−5 for all z∈S.Three functions are used(45)f1x,y=5,x,y∈D,f2x,y=x2+y2,x,y∈D,f3x,y,k=expikr,r=x2+y2,x,y∈D,where i=−1. We consider two domains (supports of the source functions) D1 and D2 (see Figure 1), the unit square with vertices(46)−1,−3,0,−3,0,−2,−1,−2,and an L-shaped domain with vertices(47)2,2,3,2,3,2.25,2.25,2.25,2.25,3,2,2.We first consider the constant functionf1=5. In Figure 2, we plot the contours of the indicator function Iz=gzα in the sampling domain. It is clear that the values of Iz is small in the neighborhood of the source. The boundary of the exact support is the solid line. The circle centered at the minimum value of Iz (dotted line) is the output of the two-step ESM, which correctly provides the location and rough size of the source for both wave numbers k1=1 and k2=6.Figure 2
Reconstructions forf1=5. (a) D1, k=1. (b) D1, k=6. (c) D2, k=1. (d) D2, k=6.
(a)
(b)
(c)
(d)Figure3 is the results for f2=x2+y2; that is, the source is a function of position.Figure 3
Reconstructions forf2=x2+y2. (a) D1, k=1. (b) D1, k=6. (c) D2, k=1. (d) D2, k=6.
(a)
(b)
(c)
(d)Finally, in Figure4, the results for f3x,y,k=expikr are shown. The reconstructions again correctly provide the location and rough size of the source, and the dependence of f3 on k does not change the results significantly.Figure 4
Reconstructions forf3x,y=expikr. (a) D1, k=1. (b) D1, k=6. (c) D2, k=1. (d) D2, k=6.
(a)
(b)
(c)
(d)


## 5. Conclusions

A two-step ESM is proposed for the inverse source problem. The ESM was originally developed for the inverse scattering problem. The new setup of the far-field equation makes it possible to treat the inverse source problem. The behavior of the solutions is theoretically justified. Numerical examples show that the two-step ESM can be used to find the location and rough size of the support of the acoustic source effectively.



---

*Source: 6434607-2020-02-07.xml*