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
dyzhang@jlu.edu.cn ykguo@hit.edu.cn li.jz@sustc.edu.cn hongyu.liuip@gmail.com
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

(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

(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

(2.3)
(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

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

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

(2.5)

where

(2.6)

and

(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

(2.8)

where .

Proof. Denote and , then we have

(2.9)

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

(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

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

(2.11)
(2.12)

and

(2.13)

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

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:

Lemma 2.2

For , we have

(2.14)
(2.15)
(2.16)

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

(2.17)
(2.18)

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

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

(2.19)

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

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

which leads to the expansion (2.18).

Let and . Then we have

which, together with Lemma 2.1, yields

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

and

where . Hence, we obtain

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

(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

Lemma 2.4

For , we have

(2.21)
(2.22)
(2.23)

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

(2.24)
(2.25)

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

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

(2.26)

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

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

which leads to the expansion (2.25).

Let and . Then we have

which, together with Lemma 2.3, yields

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

and

where . Hence, we obtain

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

(2.27)
(2.28)

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

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.

Algorithm DSM: Direct sampling method for recovering point sources
Step 1 Collect the Cauchy data and with the fixed wavenumber ;
Step 2 Select a sampling grid over the probe region ;
Step 3 For each sampling point , compute the values of the indicator functions ;
Step 4 According to the images of , collect all the significant local maximizers , where denotes the estimated number of sources. If several significant local maximizers are clustered in a region whose diameter is far less than , then they are treated as a single local maximizer;
Step 5 Cluster , into groups such that in each group the distances between distinct points are less than ;
Step 6 Take the average location of all the points in each group as the reconstruction of the locations of the sources;

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.

Algorithm DSM2: Two-level DSM for recovering point sources
Step 1 Collect the Cauchy data and with the fixed wavenumber ;
Step 2 Select a global coarse sampling grid over the probe region ;
Step 3 For each sampling point , compute the values of the indicator functions ;
Step 4 According to the images of , collect all the significant local maximizers , where denotes the estimated number of sources. If several significant local maximizers are clustered in a region whose diameter is far less than , then they are treated as a single local maximizer;
Step 5 For each local maximizer , select a local fine sampling grid centered at such that its diameter is less than ;
Step 6 Locate the maximizers of over for ;
Step 7 Cluster , into groups such that in each group the distances between distinct points are less than ;
Step 8 Take the average location of all the points in each group as the reconstruction of the locations of the sources;
(a)
(b)
Figure 1: An illustration of DSM2 in 2D. (a) global coarse grid (black lines). The maximizers of the indicator function over the global coarse grid are denoted by the small red points; (b) local fine grid (blue lines). In the vicinity of each maximizer, a local fine grid is used to resample the concerned region and enhance the accuracy of the maximizer location.
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

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