Locating multiple multipolar acoustic sources using the direct sampling method

# Locating multiple multipolar acoustic sources using the direct sampling method

Deyue Zhang, Yukun Guo, Jingzhi Li and Hongyu Liu School of Mathematics, Jilin University, Changchun, P. R. China
Department of Mathematics, Harbin Institute of Technology, Harbin, P. R. China
Department of Mathematics, Southern University of Science and Technology, Shenzhen, P. R. China
Department of Mathematics, Hong Kong Baptist University,Kowloon Tong, Hong Kong SAR
###### Abstract

This work is concerned with the inverse source problem of locating multiple multipolar sources from boundary measurements for the Helmholtz equation. We develop simple and effective sampling schemes for location acquisition of the sources with a single wavenumber. Our algorithms are based on some novel indicator functions whose indicating behaviors could be used to locate multiple multipolar sources. The inversion schemes are totally “direct” in the sense that only simple integral calculations are involved in evaluating the indicator functions. Rigorous mathematical justifications are provided and extensive numerical examples are presented to demonstrate the effectiveness, robustness and efficiency of the proposed methods.

, , ,

## 1 Introduction

The inverse source problem concerns with the reconstruction of unknown sources from measured scattering data away from the sources. It arises in many scientific fields and engineering applications, such as antenna synthesis [6, 16, 39], active acoustic tomography [5, 40], medical imaging [3, 4, 7, 25] and pollution source tracing [21, 26].

Over recent years, intensive attention [6, 8, 9, 16, 22, 23, 24, 36, 37, 38, 39] has been focused on the inverse source problem of determining a source in the Helmholtz equation

 Δu+k2u=Fin Ω, (1.1)

from boundary measurements and , where is the wavenumber, () is a bounded Lipschitz domain with boundary and denotes the outward unit normal to . A main difficulty of the inverse source problem with a single wavenumber is the non-uniqueness of the source due to the existence of non-radiating sources [3, 10, 11, 15, 17], and several numerical methods with multi-frequency measurements [8, 9, 24, 42] have been proposed to overcome it for the source with a compact support in the sense. However, fortunately, with a single wavenumber, the uniqueness can be obtained if a priori information on the source is available [18, 23]. Readers may refer to [18, 19, 20, 27] and the reference therein for the further stability results.

In this paper, we assume that the source is a finite combination of well separated monopoles and dipoles of the form

 F(x)=M∑j=1(λj+ηj⋅∇)δ(x−zj), (1.2)

where stands for the Dirac distribution, signifies the number of the source points, are points in , and and are respectively, scalar and vector source intensities such that and . Here, the points are assumed to be mutually distinct. The inverse source problem under concern is a location acquisition problem, which can be stated as: Find the locations of the source of the form (1.2) in the Helmholtz equation (1.1) from the boundary measurements and with a single wavenumber .

Numerical methods for determining the multipolar sources have drawn a lot of interest in the literature. For the Poisson equation, an algebraic method for recovering monopolar sources () was proposed in [22], and has been extended to the case of multipolar sources in [12, 13, 34, 35]. The algebraic method has also been developed to the 3D Helmholtz equation in [23] for monopolar sources, and then to the 3D elliptic equations in [2] for multipolar sources. In the case of 2D elliptic equations, the method has been extended in [1] for monopolar sources. We refer to [28] for a relevant paper on the reconstruction of extended sources for the 2D Helmholtz equation.

The purpose of this paper is to provide simple and effective numerical methods for determining the multipolar sources for the Helmholtz equation in two and three dimensions. We focus our attention on the non-iterative sampling-type methods, and we shall develop some novel direct sampling schemes in this paper. The main technicality lies in the decaying property of oscillatory integrals, which is employed to construct indicator functions at any sampling point such that the proposed indicator functions attain the local maximum at in an open neighborhood of . We would like to emphasize that the proposed methods follow the one-shot ideas [29, 30, 32, 33] and possess the following attractive features: (a) the indicator functions can be formulated in close form via the boundary data; (b) the imaging schemes are explicit since the imaging indicator functions do not rely on any matrix inversion or forward solution process; (c) our methods are easy to implement with computational efficiency due to the fact that only cheap integrations are involved in the indicator functions; (d) the locating schemes perform robust with respect to large noise level.

The outline of this paper is as follows. In Section 2, we will present the rationale behind the new method for the 2D and 3D case, respectively. In each case, rigorous mathematical justifications of the proposed sampling schemes are provided in detail. The direct sampling methods (DSM) and its two-level version are proposed at the end. In Section 3, uniqueness and stability of the proposed methods are established and analyzed in terms of measurement error, respectively. Numerical examples are demonstrated in Section 4 to verify the performance of the proposed reconstruction schemes. Finally, we conclude the work with some remarks and discussions of future topics.

## 2 Direct sampling method

In this section, we will present our sampling schemes and give rigorous mathematical justifications. First, we make the following assumptions

 L:=min1≤j,j′≤Mj≠j′dist(zj,zj′)≫2πk, (2.3) |λj||ηj′|∼kfor λj≠0 and ηj′≠0,∀j,j′=1,⋯,M.  REMOVED!! (2.4)

Condition (2.3) means that the sources are sparsely distributed, i.e., they are well separated with respect to the wavelength . Condition (2.4) indicates that for larger wavenumber , the magnitudes of the data produced by monopoles and dipoles should be of the same order. Otherwise, for example, if , then the information from the monopoles would be annihilated by the data due to the dipoles, such as

 ∣∣∣eik|x|4π|x|∣∣∣=14π|x|and∣∣∣∇eik|x|4π|x|∣∣∣=k4π|x|+O(|x|−2);

particularly, with the noisy data, one can not expect a reasonable numerical result on the monopoles.

Let be the unit sphere in and be represented by

 d={(d1,d2),N=2,(d1,d2,d3),N=3,

For uniformity of exposition, we introduce a constant throughout this paper. For any sampling point , we define indicator functions

 IN,ℓ(z):=aN,ℓ2N−1π∫SN−1R(d)dℓe−ikd⋅zds(d),ℓ=0,1,⋯,N, (2.5)

where

 aN,ℓ:=⎧⎨⎩1,ℓ=0,Nik,ℓ=1,⋯,N, (2.6)

and

 R(d):=∫Γ(eikx⋅d∂νu(x)−u(x)∂νeikx⋅d)ds(x). (2.7)

These indicator functions will be used to find the locations of the sources.

### 2.1 Two dimensional case

To see the characteristics of the indicator functions in two dimensions, we need to establish two crucial lemmas.

###### Lemma 2.1

Let , . Then

 ∣∣∣∫S1dpdqeikd⋅zds(d)∣∣∣=O((k|z|)−1/2),k|z|→∞, (2.8)

where .

Proof. Denote and , then we have

 ∫S1eikd⋅zds(d)=∫2π0eik|z|cos(θ−α)dθ. (2.9)

Introduce the Jacobi-Anger expansion (see [41, Section 2.22] and [14, p.75])

 eiρcosϕ=+∞∑n=−∞inJn(ρ)einϕ, (2.10)

where is the Bessel function of the first kind of order . By taking and in (2.10), then integrating over with respect to and using (2.9), we obtain

 ∫S1eikd⋅zds(d)=2πJ0(k|z|).

Similarly, from , , , and (2.10), we derive

 ∫S1d1eikd⋅zds(d) = ∫2π0eik|z|cos(θ−α)cosθdθ (2.11) = +∞∑n=−∞in2Jn(k|z|)∫2π0ein(θ−α)(eiθ+e−iθ)dθ = i−12J−1(k|z|)eiα2π+i2J1(k|z|)e−iα2π = 2πicosα J1(k|z|),
 ∫S1d2eikd⋅zds(d)=2πisinα J1(k|z|),
 ∫S1d21eikd⋅zds(d)=πJ0(k|z|)−πcos2α J2(k|z|), (2.12)
 ∫S1d22eikd⋅zds(d)=πJ0(k|z|)+πcos2α J2(k|z|),

and

 ∫S1d1d2eikd⋅zds(d)=−πsin2α J2(k|z|). (2.13)

Hence, the estimate (2.8) follows from the asymptotic behavior (see [14, p.74])

 H(1)n(t)=√2πtei(t−nπ2−π4){1+O(t−1)},t→∞,

and , where denotes the Hankel function of the first kind of order .

The following lemma follows from the definition of the Bessel functions by truncating the infinite series:

 J0(t) = ∞∑p=0(−1)p(p!)2(t2)2p, Jn(t) = ∞∑p=0(−1)pp!(n+p)!(t2)n+2p = tn2nn!⎛⎜ ⎜⎝1+∞∑p=1(−1)p(t2)2pp!(n+1)⋯(n+p)⎞⎟ ⎟⎠,n≥1.
###### Lemma 2.2

For , we have

 0

Now, we present the indicating behaviors of , which play a vital role in our schemes.

###### Theorem 2.1

Let source be of the form (1.2) with satisfying and , the assumptions (2.3), (2.4) hold, and be described as in (2.5). Then, we have the following asymptotic expansions

 I2,0(zj)=λj+O((kL)−1/2), (2.17) I2,ℓ(zj)=ηj,ℓ+O((kL)−1/2),ℓ=1,2. (2.18)

Moreover, there exists an open neighborhood of , , such that

 |I2,0(z)|≤|λj|+O((kL)−1/2),z∈neigh(zj)forλj≠0, |I2,ℓ(z)|≤|ηj,ℓ|+O((kL)−1/2),z∈neigh(zj)forηj,ℓ≠0,ℓ=1,2,

where the equalities hold only at . That is, is a local maximizer of for and for in .

Proof. Without loss of generality, we only consider the indicating behavior of . First, by multiplying equation (1.1) by for , then integrating over , using Green’s formula and (2.7) we obtain

 M∑j=1λjeikd⋅zj−ikM∑j=1(ηj⋅d)eikd⋅zj=R(d),d∈S1. (2.19)

Next, we multiply equation (2.19) by for sampling point , integrate over , and then derive

 I2,1(z) = 1πM∑j=1∫S1(ηj,1d1d1+ηj,2d1d2)eikd⋅(zj−z)ds(d) +ikπM∑j=1λj∫S1d1eikd⋅(zj−z)ds(d).

Further, by using the assumptions (2.3), (2.4), and Lemma 2.1, we have for ,

 I2,1(zj′) = +ikπM∑j=1j≠j′λj∫S1d1eikd⋅(zj−zj′)ds(d) = ηj′,1+O((kL)−1/2),

which leads to the expansion (2.18).

Let and . Then we have

 |zj−z|≥L−ρ≫2πk,z∈Bρ(zj′), j≠j′,

which, together with Lemma 2.1, yields

 I2,1(z) = 1π∫S1(ηj′,1d1d1+ηj′,2d1d2)eikd⋅(zj′−z)ds(d) +ikπλj′∫S1d1eikd⋅(zj′−z)ds(d)+O((kL)−1/2),z∈Bρ(zj′).

By using (2.11)–(2.16), we deduce that for and ,

 1π∣∣∣∫S1d1d1eikd⋅(zj′−z)ds(d)∣∣∣ ≤ J0(k|zj′−z|)+J2(k|zj′−z|) < 1−τ28+τ464 < 1,
 1π∣∣∣∫S1d1d2eikd⋅(zj′−z)ds(d)∣∣∣≤J2(k|zj′−z|)<τ28

and

 1kπ∣∣∣∫S1d1eikd⋅(zj′−z)ds(d)∣∣∣≤2kJ1(k|zj′−z|)<τk,

where . Hence, we obtain

 |I2,1(z)|<|ηj′,1|+O((kL)−1/2)for z∈Bρ(zj′)∖{zj′}.

This completes the proof of the theorem.

### 2.2 Three dimension case

In this subsection, we shall analyze the properties of the indicator functions in 3D. Before showing the indicating behaviors, we establish two crucial lemmas.

###### Lemma 2.3

Let , . Then

 ∣∣∣∫S2dpdqeikd⋅zds(d)∣∣∣=O((k|z|)−1),k|z|→∞, (2.20)

where .

The proof is technical and lengthy, therefore we deter it till Appendix A.

The following lemma is a 3D version of Lemma 2.2 and follows from the definition of the spherical Bessel functions

 jn(t)=∞∑p=0(−1)ptn+2p2pp!1⋅3⋅⋅⋅(2n+2p+1).
###### Lemma 2.4

For , we have

 0

Now, the indicating behaviors of are presented in the following theorem.

###### Theorem 2.2

Let source be of the form (1.2) with satisfying and , the assumptions (2.3), (2.4) hold, and indicator functions be described as in (2.5). Then, we have the following asymptotic expansions

 I3,0(zj)=λj+O((kL)−1), (2.24) (2.25)

Moreover, there exists an open neighborhood of , , such that

 |I3,0(z)|≤|λj|+O((kL)−1),z∈neigh(zj)forλj≠0, |I3,ℓ(z)|≤|ηj,ℓ|+O((kL)−1),z∈neigh(zj)forηj,ℓ≠0,ℓ=1,2,3,

where the equalities hold only at . That is, is a local maximizer of for and for in .

Proof. Without loss of generality, we only consider the indicating behavior of . First, by multiplying equation (1.1) by for , then integrating over , using Green’s formula and (2.7) we obtain

 M∑j=1λjeikd⋅zj−ikM∑j=1(ηj⋅d)eikd⋅zj=R(d),d∈S2. (2.26)

Next, we multiply equation (2.26) by for sampling point , integrate over , and then derive

 I3,1(z) = 34πM∑j=1∫S2(ηj,1d1d1+ηj,2d1d2+ηj,3d1d3)eikd⋅(zj−z)ds(d) +3i4kπM∑j=1λj∫S2d1eikd⋅(zj−z)ds(d).

Further, by using the assumptions (2.3), (2.4), and Lemma 2.3, we have for ,

 I3,1(zj′) = ηj′,1+34πM∑j=1j≠j′∫S2(ηj,1d1d1+ηj,2d1d2+ηj,3d1d3)eikd⋅(zj−zj′)ds(d) +3i4kπM∑j=1j≠j′λj∫S2d1eikd⋅(zj−zj′)ds(d) = ηj′,1+O((kL)−1),

which leads to the expansion (2.25).

Let and . Then we have

 |zj−z|≥L−ρ≫2πk,z∈Bρ(zj′), j≠j′,

which, together with Lemma 2.3, yields

 I3,1(z) = 34π∫S2(ηj′,1d1d1+ηj′,2d1d2+ηj′,3d1d3)eikd⋅(zj′−z)ds(d) +3i4kπλj′∫S2d1eikd⋅(zj′−z)ds(d)+O((kL)−1),z∈Bρ(zj′).

By using (A.38)–(2.23), we deduce that for and ,

 34π∣∣∣∫S2d1d1eikd⋅(zj′−z)ds(d)∣∣∣ ≤ j0(k|zj′−z|) +|3cos2α−1|+sin2α2j2(k|zj′−z|) ≤ j0(k|zj′−z|)+j2(k|zj′−z|) < 1−τ210+τ4120 < 1,
 34π∣∣∣∫S2d1d2eikd⋅(zj′−z)ds(d)∣∣∣≤32j2(k|zj′−z|)<τ210,
 34π∣∣∣∫S2d1d3eikd⋅(zj′−z)ds(d)∣∣∣≤3j2(k|zj′−z|)<τ25

and

 34kπ∣∣∣∫S2d1eikd⋅(zj′−z)ds(d)∣∣∣≤3kj1(k|zj′−z|)<τk,

where . Hence, we obtain

 |I3,1(z)|<|ηj′,1|+O((kL)−1)for z∈Bρ(zj′) and z≠zj′.

This completes the proof of the theorem.

### 2.3 Sampling schemes

Now we are in the position to present the main theorem of this paper by combining Theorems 2.1 and 2.2.

###### Theorem 2.3

Let source be of the form (1.2) with satisfying and , the assumptions (2.3), (2.4) hold, and indicator functions be described as in (2.5). Then, we have the following asymptotic expansions

 IN,0(zj)=λj+O((kL)−N−12), (2.27) IN,ℓ(zj)=ηj,ℓ+O((kL)−N−12),ℓ=1,⋯,N. (2.28)

Moreover, there exists an open neighborhood of , , such that for all , it holds that

 |IN,0(z)|≤|λj|+O((kL)−N−12)forλj≠0, |IN,ℓ(z)|≤|ηj,ℓ|+O((kL)−N−12)forηj,ℓ≠0,ℓ=1,⋯,N,

where the equalities hold only at . That is, is a local maximizer of for and for in .

Motivated by Theorem 2.3, we are ready to present the first direct sampling method in , see Algorithm DSM.

Our numerical experiments showed that the sampling scheme DSM could work very well if the sampling grid is sufficiently fine. However, a uniform fine mesh is computationally expensive and it would be more efficient to use a relatively coarse sampling grid for the most part of the probe region since the point sources are assumed to be sparsely distributed. On the other hand, it is difficult to locate the local maximizers accurately if only a uniformly coarse grid is available. To guarantee the accuracy and computational efficiency simultaneously, we propose a two-level sampling strategy, whose main feature consists of a global coarse grid for rough locating and a local fine grid for fine tuning. The two-level direct sampling scheme is stated as Algorithm DSM2. To sketch the idea of DSM2, we refer to Figure 1.

###### Remark 2.1

In addition, if we know a priori that (resp. ), i.e., the target source consists of only monopoles (resp. dipoles), then according to Theorem 2.3, the above algorithms could be simplified by utilizing only indicator(s) (resp. ).

## 3 Uniqueness and stability

We begin this section with an interesting result on the uniqueness under the assumptions (2.3) and (2.4), which indicates that our sampling schemes can well distinguish the sparsely distributed sources.

###### Theorem 3.1

Let source be of the form (1.2) and the assumptions (2.3), (2.4) hold, then the source can be uniquely determined by the boundary measurements and .

Proof. Without loss of generality, we only need to consider the homogeneous boundary value problem for the 3D case. Let . Then we have , , and thus, . In particular, . Furthermore, by using (2.27) and (2.28), we obtain

 λj=O((kL)−N−12) and ηj,ℓ=O((kL)−N−12), ℓ=1,⋯,N,

which, together with assumptions (2.3) and (2.4), yields and . This completes the proof.

In the following, we are going to analyze the stability of our methods. Let the measured noisy data satisfy

 ∥uϵ−u∥L2(Γ)≤ϵ∥u∥L2(Γ),∥∂νuϵ−∂νu∥L2(Γ)≤