Numerical schemes to reconstruct three dimensional time-dependent point sources of acoustic waves

TL;DR

This paper introduces numerical schemes based on the modified method of fundamental solutions for reconstructing three-dimensional time-dependent acoustic point sources, including stationary and moving sources, demonstrating their effectiveness through numerical experiments.

Contribution

The paper proposes a novel MMFS-based approach for 3D time-dependent inverse acoustic source problems, including a simplified sampling method for moving sources.

Findings

Effective reconstruction of stationary point sources

Successful reconstruction of moving point sources

Numerical experiments validate the proposed methods

Abstract

This paper is concerned with the numerical simulation of three dimensional time-dependent inverse source problems of acoustic waves. The reconstructions of both multiple stationary point sources and a moving point source are considered. The modified method of fundamental solutions (MMFS), which expands the solution utilizing the time convolution of the Green's function and the signal function, is proposed to solve the problem. For the reconstruction of a moving point source, moreover, the MMFS is simplified as a simple sampling method at each time step. Numerical experiments are provided to show the effectiveness of the proposed methods.

Tables

Algorithm I.  MMFS1 to reconstruct stationary point sources
Step 1	Choose a convex region $\Omega$, a signal function $\lambda \left(t\right)$, an integer $M$ and the locations
	$s_j\left(j=1,\dots ,M\right)$ of the point sources. Collect the wave data $u(x_i,t_k)$ for the
	sensing points $x_i\in \partial \Omega \left(i=1,\dots ,N_x\right)$ and the discrete time steps $t_k\in [0,T]$
	$\left(k=1,\dots ,N_T\right)$, where $T$ is a chosen terminal time.
Step 2	Choose a sampling region $D\subset \Omega$ such that $s_j\in D$ and $D\cap \partial \Omega =\varnothing$. Select a grid
	of sampling points $z_l\left(l=1,\dots ,N_z\right)$ in $D$. Compute $c\left(z_l\right)$ from
	$\sum _{l=1}^{N_z}c\left(z_l\right)\left(G\ast \lambda \right)(x_i,t_k;z_l)=u(x_i,t_k),i=1,\dots ,N_x,k=1,\dots ,N_T$
	using the conjugate gradient method.
Step 3	Mesh $c\left(z_l\right)$ on the sampling grid. The locations of the point sources are given by
	the locations of $z_l$ for which $c\left(z_l\right)$ are local maximum values.

Algorithm II.  The simplified scheme to reconstruct a moving point source
Step 1	Choose a convex region $\Omega$, a signal function $\lambda \left(t\right)$ and the trajectory $s\left(t\right),t\in [0,T]$
	of the moving point source. Collect the wave data $u(x_i,t_k)$ on the sensing points
	$x_i\in \partial \Omega \left(i=1,\dots ,N_x\right)$ and the discrete time steps $t_k\left(k=1,\dots ,N_T\right)$.
Step 2	Choose a sampling region $D\subset \Omega$ such that $s\left(t\right)\subset D$ and $D\cap \partial \Omega =\varnothing$. Select a grid
	of sampling points $z_l\left(l=1,\dots ,N_z\right)$ in $D$. For each time step $t_k$, compute
	$I(z_l,t_k)={\left(\sum _{i=1}^{N_x}{\left(\left(G\ast \lambda \right)(x_i,t_k;z_l)-u(x_i,t_k)\right)}^2\right)}^{-1/2}.$
Step 3	For each time step $t_k$, the location $s\left(t_k\right)$ of the point source is approximated by
	the location of $z_l$ for which $I(z_l,t_k)$ is the global maximum value.

Table 3. Table 1: Reconstruction of the locations and intensities of the point sources.

No.	The point sources		The reconstructions
	Location	Intensity	Location	Intensity
1	$(0,1,0)$	2	$(0,0.91,0)$	2.01
2	$(-1,-1,0)$	3	$(-1.09,-1.09,0)$	5.17
3	$(-1,-1.2,0)$	2	Null	Null
4	$(1,0.5,0)$	4	$(0.91,0.55,0)$	4.31
5	$(-0.5,0.5,0)$	3	$(-0.55,0.55,0)$	3.42
6	$(1.5,-1,0)$	3	$(1.45,-0.91,0)$	3.13

Equations

$$G(x,t;s)=\frac{\delta (t-c^{-1}|x-s|)}{4\pi |x-s|},$$

$$G(x,t;s)\ast \lambda (t)=\frac{\lambda (t-c^{-1}|x-s|)}{4\pi |x-s|},$$

$$c^{-2}{\partial }_{tt}u(x,t)-\Delta u(x,t)=\lambda (t)\sum _{j=1}^M{a}_j\delta (x-s_j),x\in {\mathbb{R}}^3,t\in \mathbb{R},$$

$$u(\cdot ,0)={\partial }_{t}u(\cdot ,0)=0\text{in}{\mathbb{R}}^3$$

$$u(x,t),x\in \partial \Omega ,t\in \mathbb{R}.$$

$$d:=\underset{x\in \partial \Omega ,s\in S}{\min }|x-s|=\left|x_0-s_k\right|,$$

$$d_i^{\prime }:=\underset{s\in S\backslash \left\{s_k\right\}}{\min }\left|x_0-s\right|>d.$$

$$B_{\varsigma }(x_0):=\left\{x:\left|x-x_0\right|<\varsigma \right\}\subset V\subset \overline{\Omega },$$

$$f_i(x,t)={\lambda }^{(i)}(t)\sum _{j=1}^{M_i}{a}_j^{(i)}\delta \left(x-s_j^{(i)}\right),i=1,2$$

$$u_1=u_2\text{on}\partial \Omega \times \mathbb{R}.$$

$$c^{-2}{\partial }_{tt}w-\Delta w=f_1-f_2$$

$$w=0$$

$$w_0(x,t):=\sum _{j=1}^{M_1}{a}_j^{(1)}G\left(x,t;s_j^{(1)}\right)\ast {\lambda }^{(1)}(t)-\sum _{j=1}^{M_2}{a}_j^{(2)}G\left(x,t;s_j^{(2)}\right)\ast {\lambda }^{(2)}(t)$$

$$S^{(i)}:=\bigcup _{j=1}^{M_i}\left\{s_j^{(i)}\right\},i=1,2$$

$$d_i:=\underset{x\in \partial \Omega ,s\in S^{(i)}}{\min }|x-s|=\left|x_0^{(i)}-s_{k_i}^{(i)}\right|,i=1,2,$$

$$d_{i}^{\prime}:=\min\limits_{s\in S^{(i)}\big{\backslash}\left\{s_{k_{i}}^{(i)}\right\}}\left|x_{0}^{(i)}-s\right|>d_{i},\quad\quad i=1,2.$$

$$\{\begin{array}{ll}{\lambda }^{(i)}(t)=0,& t\le t_i,\\ {\lambda }^{(i)}(t)\ne 0,& t_i<t<t_i+\tau ,\end{array}i=1,2$$

$$w_0\left(x_0^{(1)},t\right)=a_{k_1}^{(1)}\frac{{\lambda }^{(1)}(t-c^{-1}{d}_1)}{4\pi d_1},t\in (T_1,T_2),$$

$$t_1+c^{-1}{d}_1=t_2+c^{-1}{d}_2.$$

$$d_3:=\underset{s\in S^{(2)}}{\min }\left|x_0^{(1)}-s\right|.$$

$$t_1+c^{-1}{d}_1=t_2+c^{-1}{d}_3.$$

$$\left|x_0^{(1)}-s_{k_2^{\prime \prime }}^{(2)}\right|=d_3=d_2=\underset{x\in \partial \Omega ,s\in S^{(2)}}{\min }|x-s|.$$

$$d_{3}^{\prime}:=\min\limits_{s\in S^{(2)}\big{\backslash}\left\{s_{k_{2}^{\prime\prime}}^{(2)}\right\}}\left|x_{0}^{(1)}-s\right|.$$

$$\left|x^{\ast }-s_{k_2^{\prime \prime }}^{(2)}\right|-\left|x^{\ast }-s_{k_1}^{(1)}\right|<d_3-d_1.$$

$$w_0(x^{\ast },t)=-a_{k_2^{\prime \prime }}^{(2)}\frac{{\lambda }^{(2)}\left(t-c^{-1}\left|x^{\ast }-s_{k_2^{\prime \prime }}^{(2)}\right|\right)}{4\pi \left|x^{\ast }-s_{k_2^{\prime \prime }}^{(2)}\right|},t\in (T_3,T_4),$$

$$\zeta (t)=a_{k_1}^{(1)}{\lambda }^{(1)}(t)-a_{k_2^{\prime \prime }}^{(2)}{\lambda }^{(2)}(t),$$

$$\{\begin{array}{ll}\zeta (t)=0,& t\le t_3,\\ \zeta (t)\ne 0,& t_3<t<t_3+{\tau }^{\prime }\end{array}$$

$$w_0(x,t)=\sum _{\begin{array}{c}j=1,...,M_1\\ j\ne k_1\end{array}}{a}_j^{(1)}G\left(x,t;s_j^{(1)}\right)\ast {\lambda }^{(1)}(t)-\sum _{\begin{array}{c}j=1,...,M_2\\ j\ne k_2^{\prime \prime }\end{array}}{a}_j^{(2)}G\left(x,t;s_j^{(2)}\right)\ast {\lambda }^{(2)}(t).$$

$$w_0(x,t)=\sum _{j=1}^{M_1-M_2}{a}_j^{(1^{\prime })}G\left(x,t;s_j^{(1^{\prime })}\right)\ast {\lambda }^{(1)}(t),$$

$$u(x,t)=\sum _{l=1}^{N_z}c(z_l)G(x,t;z_l)\ast \lambda (t),$$

$$c^{-2}{\partial }_{tt}u(x,t)-\Delta u(x,t)=\lambda (t)\sum _{j=1}^M{a}_j\delta (x-s_j),x\in {\mathbb{R}}^3,t\in \mathbb{R},$$

Full text

Keywords: time-dependent, inverse source problem, wave equation, modified method of fundamental solutions, sampling method

1 Introduction

The inverse problems for partial differential equations appear in various fields of science and engineering, and have been extensively studied in the past decades [18, 22, 24]. Among them, the inverse source problem, especially the identification of moving sources, has a wide range of applications such as under water sonar [27, 28], sound simulation and sound source localization [19, 29].

For the reconstruction of stationary sources, the inverse source problems with sources $\delta(t)g(x)$ that are delta-like in time and of limited oscillation in space, and sources $q(t)\delta(\partial G)$ that are oscillation in time and delta-like on the boundary of a region are considered in [13] and [34], respectively. Uniqueness analysis related to the Helmholtz equation with phaseless data is shown in [38, 39]. The stability analysis and identification of multiple point sources for the time-harmonic case are considered in [3, 4]. The conditional stability estimate of the wave equation on a line related to the inverse source problem is provided by [11, 12]. Multi-frequency inverse source problems are analyzed in [7, 8, 26, 37]. Analysis of random sources can be seen in [5, 25]. Time-dependent inverse source problems in elastodynamics are analyzed in [6].

For the reconstruction of a moving point source, direct identifications of the moving point source are studied in [30, 35]. Analysis of the moving point source when the velocity of the source is comparable to the speed of wave propagation can be seen in [15]. Matched-filter imaging method and correlation-based imaging for small fast moving debris with constant velocity are analyzed in [14]. A gesture-based input technique with the electromagnetic wave is analyzed in [21].

The method of fundamental solutions (MFS) is a meshless method which expands the solution utilizing the fundamental solution [2, 10, 32, 36]. The property of the fundamental solution, or the Green’s function, is the theoretical basis of the MFS. However, the Green’s function of the d’Alembert operator $c^{-2}\partial_{tt}-\Delta$ is

[TABLE]

where $c>0$ denotes the sound speed of the homogeneous background medium, $\partial_{tt}u=\frac{\partial^{2}u}{\partial t^{2}}$, $\Delta$ is the Laplacian in $\mathbb{R}^{3}$, and $\delta$ is the Dirac delta distribution. Since the Green’s function involves the Dirac delta distribution, the MFS is no longer feasible to solve the three dimensional wave equation. Unable to be applied directly, the Green’s function of the d’Alembert operator usually appears in the time convolution

[TABLE]

where $\lambda(t)$ is a signal function. One of the most popular application of $G(x,t;s)*\lambda(t)$ is that in the boundary integral equation method, which is a commonly used method [9, 17, 31, 33]. Therefore, instead of the Green’s function, new bases $G(x,t;s)*\lambda(t)$ are employed in the modified method of fundamental solutions (MMFS) proposed in this paper. Moreover, the MMFS can be simplified to a simple sampling method at each time step to reconstruct a moving point source, in which the sampling method is a well-known method in the numerical computation of inverse problems [16, 20, 22, 23, 40].

The time convolution of the Green’s function and the signal function is an invaluable tool for the analysis of the time domain scattering problems. Therefore, the MMFS has important significance in the theory of the time domain analysis. Moreover, the proposed methods are feasible to reconstruct both multiple stationary point sources and a moving point source. The numerical implementations of the proposed methods are simple, and extensive experiments are provided to show the effectiveness of the methods.

The outline of this paper is as follows. In Section 2, the inverse source problem with multiple stationary point sources is considered. The uniqueness result is provided and the MMFS is proposed. In Section 3, the MMFS is applied to the inverse source problem with a moving point source. Moreover, the method is simplified as a simple sampling method at each discrete time. In Section 4, numerical experiments are provided to show the effectiveness of the proposed methods. The conclusion remarks are given in Section 5.

2 Reconstruction of stationary point sources

Denote by $\Omega\subset\mathbb{R}^{3}$ a bounded convex open region. Consider the wave equation

[TABLE]

where $M\in\mathbb{N}^{*}$ is a positive integer, $s_{j}\in\Omega$ are stationary source points, and $a_{j}>0$ are the intensities of the sources.

The source points $s_{j}$ are assumed to be mutually distinct. The signal function $\lambda(t)$ is assumed to be causal, which means $\lambda(t)=0$ for $t<0$. Thus the source term $f(t)=0$ for $t<0$, and the initial condition

[TABLE]

is a direct conclusion of the causality.

The inverse source problem (P1) under consideration is: Determine the locations and intensities of the stationary point sources in (1) from the measurement data

[TABLE]

The following lemma is needed to prove the uniqueness of the solution to the inverse source problem (P1).

Lemma 2.1**.**

Let $S:=\bigcup\limits_{j=1}\limits^{M}\left\{s_{j}\right\}$ be a set of points in a bounded convex open region $\Omega\subset\mathbb{R}^{3}$, where $M\in\mathbb{N}^{*}$. Assume that

[TABLE]

where $x_{0}\in\partial\Omega$ and $1\leq k\leq M$. Then

[TABLE]

Proof.

It is obvious that $d^{\prime}\geqslant d$. By reduction to absurdity, assume that $d^{\prime}=\left|x_{0}-s_{k^{\prime}}\right|=d$ for some $s_{k^{\prime}}\in S\backslash\left\{s_{k}\right\}$. As is shown in Figure 1, construct a rectangle $s_{k}s_{k^{\prime}}AB$ with $s_{k^{\prime}}A=s_{k}B=d$, such that $x_{0}\in s_{k}s_{k^{\prime}}AB$. Then $d=\min\limits_{x\in\partial\Omega,s\in S}|x-s|$ implies that $A,B\in\overline{\Omega}$. Thus $s_{k}s_{k^{\prime}}AB\subset\overline{\Omega}$ since $\overline{\Omega}$ is a closed convex region. Spinning $s_{k}s_{k^{\prime}}AB$ around the segment $s_{k}s_{k^{\prime}}$, we get a cylinder $V\subset\overline{\Omega}$. Then there exists a $\varsigma>0$ such that

[TABLE]

which is a contradiction to $x_{0}\in\partial\Omega$. ∎

Then we have the following uniqueness result.

Theorem 2.2**.**

Assume that $\Omega\subset\mathbb{R}^{3}$ is a bounded convex open region. Let

[TABLE]

be two source terms with $M_{i}\in\mathbb{N}^{*}$, $s_{j}^{(i)}\in\Omega$, $a_{j}^{(i)}>0$ and $\lambda^{(i)}(t)\in C(\mathbb{R})$, such that the corresponding solutions to (1) for $f_{1}$ and $f_{2}$ are $u_{1}$ and $u_{2}$, respectively. Assume that $\lambda^{(i)}(t)$ are nontrivial causal functions and

[TABLE]

Then $M_{1}=M_{2}=M$, $s_{j}^{(1)}=s_{\pi(j)}^{(2)}$ and $a_{j}^{(1)}\lambda^{(1)}(t)=a_{\pi(j)}^{(2)}\lambda^{(2)}(t)$, $j=1,2,\ldots,M$ for some permutation $\pi(j)$ of $1,2,\ldots,M$.

Proof.

Denote $w:=u_{1}-u_{2}$. Then

[TABLE]

Note that

[TABLE]

is the unique causal solution (refer to Section 1.4 of [31] for the uniqueness) of the wave equation (4).

For the convenience of the expression, the rest of the proof is divided into five parts.

(i) As is shown in Figure 2, denote

[TABLE]

and

[TABLE]

in which $x_{0}^{(i)}\in\partial\Omega$ and $1\leq k_{i}\leq M_{i}$. Notice that there may be several sets of points $x_{0}^{(i)}$ and $s_{k_{i}}^{(i)}$ which satisfy (6). Nevertheless, the choice of the points would not affect the following proof.

Moreover, Lemma 2.1 implies that

[TABLE]

(ii) Next, we are going to take into consideration of the causality. Since $\lambda^{(i)}(t)\in C(\mathbb{R})$ are nontrivial causal functions, we have

[TABLE]

for some $t_{i}>0$ and $\tau>0$. Denote $t_{i}$ as the “starting time” of the signal $\lambda^{(i)}(t)$.

Assume that $t_{1}+c^{-1}d_{1}<t_{2}+c^{-1}d_{2}$. Then

[TABLE]

where $T_{1}=t_{1}+c^{-1}d_{1}$ and $T_{2}=\min\{t_{1}+c^{-1}d_{1}^{\prime},t_{2}+c^{-1}d_{2}\}$. Note that $\lambda^{(1)}(t-c^{-1}d_{1})$ is nontrivial for $t\in(T_{1},T_{2})$. Then (5) implies $a_{k_{1}}^{(1)}=0$, which is a contradiction to $a_{j}^{(i)}>0$. Using the reduction to absurdity, we have $t_{1}+c^{-1}d_{1}\geq t_{2}+c^{-1}d_{2}$. Similarly, we can prove that $t_{2}+c^{-1}d_{2}\geq t_{1}+c^{-1}d_{1}$. Thus

[TABLE]

(iii) Denote

[TABLE]

Apparently $d_{3}\geq d_{2}$. Assume that $d_{3}>d_{2}$. Then (7) implies $t_{1}+c^{-1}d_{1}<t_{2}+c^{-1}d_{3}$. A similar discussion as that in (ii) leads to a contradiction, which means $d_{3}=d_{2}$. Then (7) implies

[TABLE]

Moreover, there is a point $s_{k_{2}^{\prime\prime}}^{(2)}\in S^{(2)}$ such that

[TABLE]

Therefore, we can reselect a new set of points $x_{0}^{(2)}$ and $s_{k_{2}}^{(2)}$ satisfying (6) with $x_{0}^{(2)}=x_{0}^{(1)}$ and $s_{k_{2}}^{(2)}=s_{k_{2}^{\prime\prime}}^{(2)}$.

Denote

[TABLE]

Then Lemma 2.1 implies $d_{3}^{\prime}>d_{3}$. Assuming that $t_{1}>t_{2}$, (8) implies $d_{1}<d_{3}$. For a point $x^{*}\in\left\{x\in\partial\Omega:0<\left|x-x_{0}^{(1)}\right|<\frac{1}{2}\min\{d_{1}^{\prime}-d_{1},d_{3}^{\prime}-d_{3}\}\right\}$, we have

[TABLE]

Then (8) implies $t_{2}+c^{-1}\left|x^{*}-s_{k_{2}^{\prime\prime}}^{(2)}\right|<t_{1}+c^{-1}\left|x^{*}-s_{k_{1}}^{(1)}\right|$ and

[TABLE]

where $T_{3}=t_{2}+c^{-1}\left|x^{*}-s_{k_{2}^{\prime\prime}}^{(2)}\right|$ and $T_{4}=\min\left\{t_{1}+c^{-1}\left|x^{*}-s_{k_{1}}^{(1)}\right|,t_{2}+\frac{1}{2}c^{-1}(d_{3}+d_{3}^{\prime})\right\}$. A similar discussion as that in (ii) implies a contradiction, which means $t_{1}\leq t_{2}$. Similarly we have $t_{2}\leq t_{1}$. Then $t_{1}=t_{2}$ and $d_{1}=d_{3}$. Referring to the proof of Lemma 2.1, we finally get $s_{k_{1}}^{(1)}=s_{k_{2}^{\prime\prime}}^{(2)}$.

(iv) Define

[TABLE]

which is the real signal function of the source point $s_{k_{1}}^{(1)}=s_{k_{2}^{\prime\prime}}^{(2)}$. Then we assert that $\zeta(t)=0$. Otherwise, we can prove that

[TABLE]

for some $t_{3}\geq t_{2}+c^{-1}\min\{d_{1}^{\prime}-d_{1},d_{3}^{\prime}-d_{3}\}>t_{2}=t_{1}$ and $\tau^{\prime}>0$. Notice that the set of source points $S:=S^{(1)}\cup S^{(2)}$ can be divided into two categories: $S^{(3)}$ such that the “starting time” of the signal function of $s\in S^{(3)}$ is $t_{3}$ and $S^{(4)}=S\setminus S^{(3)}$ such that the “starting time” of the signal function of $s\in S^{(4)}$ is $t_{1}$. If $S^{(4)}=\varnothing$, it is easy to get a contradiction. For $S^{(4)}\neq\varnothing$, a similar discussion as that in (i)-(iii) implies $t_{3}=t_{1}$, which is a contradiction to $t_{3}>t_{1}$.

Since we have proved $s_{k_{1}}^{(1)}=s_{k_{2}^{\prime\prime}}^{(2)}$ and $a_{k_{1}}^{(1)}\lambda^{(1)}(t)=a_{k_{2}^{\prime\prime}}^{(2)}\lambda^{(2)}(t)$, the wave field can be rewritten as

[TABLE]

(v) Assume that $M_{1}\neq M_{2}$, there is no harm to suppose that $M_{1}>M_{2}$. Following the procedure of (i)-(iv), we can get

[TABLE]

where $a_{j}^{(1^{\prime})}>0$. Again, it is easy to get a contradiction. Then we have $M_{1}=M_{2}=M$. Following the procedure of (i)-(iv), we can finally get $s_{j}^{(1)}=s_{\pi(j)}^{(2)}$ and $a_{j}^{(1)}\lambda^{(1)}(t)=a_{\pi(j)}^{(2)}\lambda^{(2)}(t)$, $j=1,2,\ldots,M$ for some permutation $\pi(j)$ of $1,2,\ldots,M$. ∎

The classic MFS expands the solution utilizing the Green’s function (refer to [32, 33, 36]). However, since the Green’s function of the d’Alembert operator involves the Dirac delta distribution, the MFS is no longer feasible to solve the three dimensional wave equation. Hence, consider the expansion

[TABLE]

where $N_{z}\in\mathbb{N}^{*}$, $z_{l}\in\Omega$ are the sampling points, and $c(z_{l})$ are unknown coefficients to be computed.

The expansion (9) leads to the first modified method of fundamental solutions (MMFS1). We introduce the following proposition concerning the MMFS1.

Proposition 2.3**.**

Assume that $\Omega\subset\mathbb{R}^{3}$ is a bounded convex open region. Let $u(x,t)$ be a causal wave field which solves

[TABLE]

where $M\in\mathbb{N}^{*}$, $s_{j}\in\Omega$, $a_{j}>0$, and $\lambda(t)\in C(\mathbb{R})$ is a non-trivial causal signal function. Assuming that the sampling points $z_{l}\in\Omega,\,l=1,2,\ldots,N_{z}$ and a group of corresponding constants $c(z_{l})$ satisfy

[TABLE]

where $N_{z}\in\mathbb{N}^{*}$. Denote $S_{d}:=\{s_{j}\}_{j=1}^{M}$ and $Z_{d}:=\{z_{l}\}_{l=1}^{N_{z}}$. Then $S_{d}\subset Z_{d}$. Moreover,

[TABLE]

Proof.

Notice that

[TABLE]

is the unique causal solution of the wave equation (10). Meanwhile, the wave field

[TABLE]

is a causal solution of

[TABLE]

Moreover, (11) implies

[TABLE]

That is, $f_{1}(x,t)=\lambda(t)\sum\limits_{j=1}\limits^{M}a_{j}\delta(x-s_{j})$ and $f_{2}(x,t)=\lambda(t)\sum\limits_{l=1}\limits^{N_{z}}c(z_{l})\delta(x-z_{l})$ are two solutions of the same inverse source problem. Then the conclusion follows from Theorem 2.2. ∎

Remark 2.1**.**

It is a strong hypothesis that the chosen sampling points $z_{l}\in\Omega$ and the constants $c(z_{l})$ solves (11) for $x\in\partial\Omega$, $t\in\mathbb{R}$. Nevertheless, the numerical experiments in Section 4 show the effectiveness of the MMFS1 even if the hypothesis is not satisfied.

The MMFS1 to solve the inverse source problem (P1) is shown in Algorithm I. The numerical application of Algorithm I can be seen in Section 4.

3 Reconstruction of a moving point source

In this section, the inverse source problem with a moving point source is considered. The wave equation is

[TABLE]

where $T>0$, $s:[0,T]\rightarrow\Omega$ signifies the smooth trajectory of the moving point source. Denote by

[TABLE]

the instaneous velocity of the point source. Again, $\lambda(t)$ is causal and the initial condition follows from the causality.

The inverse source problem (P2) is: Determine the trajectory $s(t)$ of the moving point source in (13) from the measurement data

[TABLE]

The MMFS1 is feasible to reconstruct stationary point sources. However, for a moving point source, the location of the point source changes over time. Thus the coefficients $c(z_{l})$ in the MMFS1 should also depend on the time variable. Therefore, consider a new expansion

[TABLE]

where $z_{l}\in\Omega$ are the sampling points and $c(\cdot;z_{l})$ are unknown functions depending on $z_{l}$. The second modified method of fundamental solutions (MMFS2) is based on (15). The algorithm of MMFS2 is similar to Algorithm I except that the new expansion (15) is employed and $c(t_{k};z_{l}),\,l=1,\ldots,N_{z}$ should be computed respectively for each time step $t_{k}$, $k=1,\ldots,N_{T}$.

Notice that there is only one point source in this case. If $|v|=0$, we have $s(t_{k})\equiv s_{0}$ for some $s_{0}\in\Omega$. Then Proposition 2.3 implies

[TABLE]

Then we expect that $G(x,t;s(t_{k}))*\lambda(t)$ is the approximation of $u(x,t_{k})$ when $|v|$ is small. On this basis, define the indicator function

[TABLE]

where $\|\cdot\|_{\partial\Omega}$ is the $L^{2}(\partial\Omega)$ norm with respect to $x$. We have the following theorem concerning the indicator function (16).

Theorem 3.1**.**

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded convex open region. Assume that $d_{\Omega}:=\sup\limits_{x,y\in\Omega}|x-y|\ll c$, $|v|\ll c$ and $\lambda,s\in C^{1}[0,T]$. Let $u_{0}(x,t)$ be the causal solution of the wave equation (13). For any fixed $t\in[0,T]$, the indicator function (16) satisfies

[TABLE]

Proof.

Note that when $|v(t)|<c$ for $t\in(0,T)$, the explicit solution to equation (13) is given by (refer to [30])

[TABLE]

where the retarded time $\tau$ satisfies $t-\tau=c^{-1}|x-s(\tau)|$.

Under the assumptions $d_{\Omega}\ll c$, $|v|\ll c$ and $\lambda,s\in C^{1}[0,T]$, we assert that

[TABLE]

is an approximation of the solution (17). The proof of a similar conclusion can be seen in [35]. Though we use an arbitrary causal signal function $\lambda(t)$ instead of the time-harmonic signal $\lambda(t)=\sin(\omega_{0}t)$ for some $\omega_{0}>0$, and the function $\lambda(t-c^{-1}|x-s(t)|)$ is occupied in $G(x,t;s(t))*\lambda(t)$ in this paper instead of $\lambda(t)$, a similar discussion implies

[TABLE]

with some $0<\varepsilon(t)\ll 1$.

Then the smoothness of the function $G(x,t;z)*\lambda(t)$ with respect to $z$ implies the conclusion. ∎

Then the MMFS2 is in fact equivalent to a simple sampling method. The simplified scheme is shown in Algorithm II. The numerical implement of Algorithm II is shown in Section 4.

4 Numerical examples

In this section, we consider the numerical implementation of the proposed algorithms. The radiated field is collected for $t\in[0,T]$, where $T$ is the terminal time. The time discretization is

[TABLE]

where $N_{T}\in\mathbb{N}^{*}$. Random noises are added to the data with

[TABLE]

where $\epsilon>0$ is the noise level and $r$ are uniformly distributed random numbers in $[-1,1]$.

In all the experiments, the signal function $\lambda(t)$ is chosen as

[TABLE]

The signal function $\lambda(t)$ and its Fourier spectrum can be seen in Figure 3.

4.1 Reconstruction of multiple static point sources

Algorithm I is employed for the reconstruction of multiple static point sources. The synthetic data $u(x,t)$, $x\in\partial\Omega$, $t\in[0,T]$ is given by the analytic solution (12). We choose $c=1$, $T=15$ and $N_{T}=64$ in this subsection.

Example 1. In this example, the reconstruction of the stationary point sources located at $(0,-1,1)$, $(0,1,-1)$, $(1,-1,0)$ and $(1,0,-1)$ with the same intensity is considered. The sensing points are chosen as

[TABLE]

with $\varphi_{i}=\frac{2i-1}{16}\pi$, $i=1,2,\ldots,8$ and $\theta_{j}=\frac{j}{4}\pi$, $j=0,1,\ldots,7$. The sampling points are $21\times 21\times 21$ uniform discrete points in $[-2,2]\times[-2,2]\times[-2,2]$. The reconstructions are shown in Figure 4.

Remark 4.1**.**

To facilitate the 3D visualization, we add 2D projections in some of the 3D figures in this paper.

Example 2. We investigate the reconstruction of point sources with different intensities in this example. The source points are chosen as $(0,1,0)$, $(-1,-1,0)$, $(-1,-1.2,0)$, $(1,0.5,0)$, $(-0.5,0.5,0)$ and $(1.5,-1,0)$ with relative intensities $2$, $3$, $2$, $4$, $3$ and $3$, respectively. The sampling points are chosen as $23\times 23\times 23$ uniform discrete points in $[-2,2]\times[-2,2]\times[-2,2]$. The sensors are chosen as all the sensors in (18), the left half of the sensors with $i=1,2,\ldots,8$, $j=0,1,\ldots,3$ and the upper half with $i=1,2,\ldots,4$, $j=0,1,\ldots,7$, respectively in three cases. The reconstructions are shown in Figure 5.

As is shown in Figure 5, the proposed method is feasible to reconstruct point sources with different intensities. The specific data of the reconstruction is given by the following procedure:

(1) Compute $c(z_{l})$ and save the data as $C(i,j,k),\,i,j,k=1,...,23$, corresponding to the sampling points $z(i,j,k)=\left(-2+\frac{2(i-1)}{11},-2+\frac{2(j-1)}{11},-2+\frac{2(k-1)}{11}\right)$, respectively. Denote $n=1$.

(2) Find a global maximum of $C(i,j,k)$ and the corresponding maximum point $z(i_{n},j_{n},k_{n})$. The intensity of the point $z(i_{n},j_{n},k_{n})$ is given by

[TABLE]

(3) If $T(i_{n},j_{n},k_{n})<1$, end the procedure.

(4) If $T(i_{n},j_{n},k_{n})>1$, the corresponding point $z_{l}(i_{n},j_{n},k_{n})$ is regarded as a source point with the intensity $T(i_{n},j_{n},k_{n})$.

(5) Denote $C(i,j,k)=0$ for $i=i_{n}-1,...,i_{n}+1$, $j=j_{n}-1,...,j_{n}+1$ and $k=k_{n}-1,...,k_{n}+1$. Redefine $n=n+1$ and go back to step (2).

The specific data is shown in Table 1. The error of the location is mainly caused by the discretization precision of the sampling region. Since the point sources No. 2 and No. 3 are too close to each other, only one source point is reconstructed with the superimposition of the intensities.

4.2 Reconstruction of a moving point source

This subsection is concerned with the reconstruction of a moving point source. Numerical scheme based on Algorithm II is employed. The synthetic data $u(x,t)$, $x\in\partial\Omega$, $t\in[0,T]$ is given by the analytic solution (17). We choose $c=340$ and $T=2\pi$ in this subsection.

Example 3. In this example, we consider the reconstruction of arbitrary trajectory of a moving source in $\mathbb{R}^{3}$. The sensors are chosen the same as that in Example 1. The sampling points are chosen as $51\times 51\times 51$ uniform discrete points in $[-3,3]\times[-3,3]\times[-3,3]$. We choose $N_{T}=64$ in this experiment. The reconstructions $s^{\prime}(t_{k})$ of the locations $s(t_{k})$ for $k=1,2,...,N_{T}$ are considered.

The trajectories of the moving source are chosen as $s_{1}(t)=(2+0.3\cos 3t)(\cos t,\sin t,0)$ and $s_{2}(t)=2(\sin{2t},\cos{2t},\frac{t}{\pi}-1)$, respectively in two cases. The reconstructions can be seen in Figure 6 and Figure 7, respectively.

As is shown in Figure 6(b), the reconstructions $s^{\prime}_{1}(t_{k})$ is close to the trajectory $s_{1}(t)$ except for several discrete points. The error given by the Euclidean distance $|s_{1}(t_{k})-s^{\prime}_{1}(t_{k})|$ at each discrete time steps $t_{k}$ can be seen in Figure 6(c). The error becomes large when the signal intensity $\lambda(t_{k})$ is near zero. Therefore, the following modification is provided after the reconstruction:

(1) If $|\lambda(t_{1})|<10^{-4}$, redefine $s^{\prime}(t_{1})=2s^{\prime}(t_{2})-s^{\prime}(t_{3})$.

(2) If $|\lambda(t_{N_{T}})|<10^{-4}$, redefine $s^{\prime}(t_{N_{T}})=2s^{\prime}(t_{N_{T}-1})-s^{\prime}(t_{N_{T}-2})$.

(3) If $|\lambda(t_{k})|<10^{-4}$ for any $k=2,\ldots,N_{T}-1$, redefine $s^{\prime}(t_{k})=\frac{1}{2}(s^{\prime}(t_{k-1})+s^{\prime}(t_{k+1}))$.

The modified reconstruction and the corresponding error are shown in Figure 6(d) and Figure 6(e), respectively. As we can see from Figure 6(e), the error $|s(t_{k})-s^{\prime}(t_{k})|$ is small at each time steps after the modification. Therefore, similar modifications are applied to all the experiments in the rest of this subsection.

The smooth reconstruction of the trajectory is given by the post-processing of the data $s^{\prime}(t_{k})$ by a Fourier approximation. The truncated Fourier expansion of order $N\in\mathbb{N}^{*}$ is employed such that

[TABLE]

where

[TABLE]

Fourier expansion of order 5 is employed to get the smooth reconstruction in this example.

Remark 4.2**.**

An important component of the error is caused by the discretization precision of the sampling region. The error $|s(t_{n})-s^{\prime}(t_{n})|<0.2$ coincides with the $51\times 51\times 51$ uniform discretization of the sampling region $[-3,3]\times[-3,3]\times[-3,3]$.

Example 4. As an addition of Example 3, the reconstruction of a handwritten Chinese character “ai” is considered. We choose $N_{T}=128$ in this experiment. The reconstructions can be seen in Figure 8.

The smooth reconstruction in this example is also provided by the Fourier expansion. However, the Chinese character “ai” has 5 strokes and can not be reconstruct with a single smooth curve. Thus the smooth reconstruction is provided respectively for each stroke. Since the point source moves faster in the gap between two strokes, we use the following strategy to provide the smooth reconstruction:

(1) If $\max\left\{|s^{\prime}(t_{k-1})-s^{\prime}(t_{k})|,|s^{\prime}(t_{k+1})-s^{\prime}(t_{k})|\right\}>0.3$ for any $k=2,\ldots,N_{T}-1$, classify $s^{\prime}(t_{k})$ as an end point of a stroke, or a point between two strokes.

(2) Separate the strokes of the character, and provide the smooth reconstruction of each stroke using the Fourier expansion.

As is shown in Figure 8, the algorithm is feasible to reconstruct the character with noise level $\epsilon=5\%$. The smooth reconstructions by the Fourier expansion with order 5 and order 3 are shown in Figure 8(c) and Figure 8(d), respectively.

Remark 4.3**.**

The smooth reconstruction by the Fourier expansion of order 5 indeed shows more details of the reconstruction than that of order 3. However, some of the details are caused by the noises. As is shown in Figure 8(c-d), the smooth reconstruction by the Fourier expansion of order 3 is better in this example.

Example 5. In this example, we are concerned about the reconstruction of the trajectory $s_{3}(t)=2\left(\frac{t}{\pi}\sin{2t},\frac{t}{\pi}\cos{2t},\frac{t}{\pi}-1\right)$ using 4 sensors. The sensing points are chosen as $(3,3,-3)$, $(3,-3,-3)$, $(-3,3,-3)$ and $(-3,-3,-3)$. The sampling points and the time discretization are chosen the same as that in Example 3. The reconstructions can be seen in Figure 9.

The error of the reconstruction with only 4 sensors is bigger than that of Example 3. Nevertheless, the smooth reconstruction ignores most of the error and the algorithm still works well.

5 Conclusion

We have considered the numerical simulation of the time dependent inverse source problems of acoustic waves. Modified methods of fundamental solutions have been established to reconstruct both multiple stationary sources and a moving point source. Moreover, the second modified method of fundamental solutions to reconstruct a moving point source has been modified to a simple sampling method. Several numerical examples have been provided to show the effectiveness of the proposed methods.

Acknowledgements

The work of Bo Chen was supported by the NSFC (No. 11671170) and the Fundamental Research Funds for the Central Universities (Special Project for Civil Aviation University of China, No. 3122018L009). The work of Yukun Guo was supported by the NSFC (No. 11601107, 41474102 and 11671111). The work of Yao Sun was supported by the NSFC (No. 11501566).