Real-time reconstruction of moving point/dipole wave sources from boundary measurements

TL;DR

This paper presents a real-time method to reconstruct multiple moving wave sources' locations and strengths from boundary measurements, using algebraic relations and reciprocity gap functionals, with demonstrated numerical effectiveness.

Contribution

It introduces a novel real-time reconstruction approach for unknown moving wave sources based on algebraic relations and reciprocity gap functionals.

Findings

Successful numerical experiments validate the reconstruction method.

Effective real-time reconstruction of source parameters from boundary data.

Method handles multiple moving sources with unknown number and properties.

Abstract

This paper is concerned with a reconstruction method for multiple moving point/dipole wave sources. We assume that the number, locations, and magnitudes/moments of wave sources are unknown, and consider the problem to reconstruct these parameters from the measurement of the wave field on the boundary. For this problem, we derive algebraic relations between the parameters of wave sources and the reciprocity gap functionals, and propose a real-time reconstruction procedure based on these relations. We perform some numerical experiments, and show the effectiveness of our reconstruction procedure.

Equations

$$\left\{\begin{array}[]{rll}\displaystyle\frac{1}{c^{2}}\partial_{t}^{2}u(t,\boldsymbol{r})-\Delta u(t,\boldsymbol{r})&=F(t,\boldsymbol{r})&\mbox{ in }(0,T)\times\Omega,\\ u(0,\boldsymbol{r})&=0,&\mbox{ in }\Omega,\\ \partial_{t}u(0,\boldsymbol{r})&=0,&\mbox{ in }\Omega,\\ u(t,\boldsymbol{r})&=0,&\mbox{ on }(0,T)\times\Gamma,\end{array}\right.$$

$$\varphi (t,\mathbf{r})={\partial }_{\nu }u(t,\mathbf{r}),\text{mbox}on[0,T]\times \Gamma ,$$

$$F(t,\mathbf{r})=\sum _{k=1}^K{q}_k(t)\delta (\mathbf{r}-{\mathbf{p}}_k(t)),$$

$$F(t,\mathbf{r})=-\sum _{k=1}^K{\mathbf{m}}_k(t)\cdot \nabla \delta (\mathbf{r}-{\mathbf{p}}_k(t)).$$

$$⟨\delta (\cdot -\mathbf{p}),v(\cdot ){⟩}_{\Omega }=v(\mathbf{p}),\mathbf{p}\in \Omega ,$$

$$⟨\nabla \delta (\cdot -\mathbf{p}),v(\cdot ){⟩}_{\Omega }=$$

$$=$$

$$\frac{1}{c^2}{\int }_{\Omega }{\partial }_{t}u(T,\mathbf{r})v(T,\mathbf{r})dV(\mathbf{r})-\frac{1}{c^2}{\int }_{\Omega }u(T,\mathbf{r}){\partial }_{t}v(T,\mathbf{r})dV(\mathbf{r})$$

$$-{\int }_0^T{\int }_{\Gamma }{\partial }_{\nu }u(t,\mathbf{r})v(t,\mathbf{r})dS(\mathbf{r})dt$$

$$+{\int }_0^T{\int }_{\Omega }u(t,\mathbf{r})\left(\frac{1}{c^2}{\partial }_t^{2}v(t,\mathbf{r})-\Delta v(t,\mathbf{r})\right)dV(\mathbf{r})dt$$

$$=$$

$$\mathcal{F}(v)=$$

$$=$$

$$\mathcal{F}(v)=$$

$$=$$

$$u{|}_{[0,T]\times (\Omega \backslash D)}\in$$

$$\nabla u{|}_{[0,T]\times (\Omega \backslash D)}\in$$

$$\nabla {\partial }_t^{m}u{|}_{[0,T]\times (\Omega \backslash D)}\in$$

$${\partial }_{\nu }u\in H^{\ell }([0,T]\times \Gamma ).$$

$$\left\{\begin{array}[]{rll}\displaystyle\frac{1}{c^{2}}\partial_{t}^{2}u_{\mathrm{N}}-\Delta u_{\mathrm{N}}&=\displaystyle\sum_{k=1}^{K}q_{k}(t)\delta(\boldsymbol{r}-\boldsymbol{p}_{k}(t)),&\qquad\mbox{in }(0,T)\times{\mathbb{R}}^{3},\\ u_{\mathrm{N}}(0,\boldsymbol{r})&=0,&\qquad\mbox{ in }{\mathbb{R}}^{3},\\ \partial_{t}u_{\mathrm{N}}(0,\boldsymbol{r})&=0,&\qquad\mbox{ in }{\mathbb{R}}^{3}.\end{array}\right.$$

$$u_{\mathrm{N}}(t,\mathbf{r})=\sum _{k=1}^K\frac{1}{4\pi }\cdot \frac{q_k(s_k(t,\mathbf{r}))}{\left|\mathbf{r}-{\mathbf{p}}_k(s_k(t,\mathbf{r}))\right|\cdot h_k(s_k(t,\mathbf{r}),\mathbf{r})},$$

$$t=s+\frac{\left|\mathbf{r}-{\mathbf{p}}_k(s)\right|}{c},$$

$$h_k(s,\mathbf{r})=1-\frac{d_t{\mathbf{p}}_k(s)\cdot (\mathbf{r}-{\mathbf{p}}_k(s))}{c\left|\mathbf{r}-{\mathbf{p}}_k(s)\right|}.$$

$$\left\{\begin{array}[]{rll}\displaystyle\frac{1}{c^{2}}\partial_{t}^{2}u_{\mathrm{H}}-\Delta u_{\mathrm{H}}&=0,&\qquad\mbox{ in }(0,T)\times\Omega,\\ u_{\mathrm{H}}(0,\boldsymbol{r})&=0,&\qquad\mbox{ in }\Omega,\\ \partial_{t}u_{\mathrm{H}}(0,\boldsymbol{r})&=0,&\qquad\mbox{ in }\Omega,\\ u_{\mathrm{H}}&=-u_{\mathrm{N}}|_{[0,T]\times\Gamma},&\qquad\mbox{ on }[0,T]\times\Gamma.\end{array}\right.$$

$$\left\{\begin{array}[]{rll}\displaystyle\frac{1}{c^{2}}\partial_{t}^{2}u_{\mathrm{N}}-\Delta u_{\mathrm{N}}&=-\displaystyle\sum_{k=1}^{K}\boldsymbol{m}_{k}(t)\cdot\nabla\delta(\boldsymbol{r}-\boldsymbol{p}_{k}(t)),&\mbox{ in }(0,T)\times{\mathbb{R}}^{3},\\ u_{\mathrm{N}}(0,\boldsymbol{r})&=0,&\mbox{ in }{\mathbb{R}}^{3},\\ \partial_{t}u_{\mathrm{N}}(0,\boldsymbol{r})&=0,&\mbox{ in }{\mathbb{R}}^{3}.\end{array}\right.$$

$$u_{\mathrm{N}}(t,\mathbf{r})=$$

$$+\sum _{k=1}^K\frac{1}{4\pi c}\cdot {\partial }_t\left(\frac{{\mathbf{m}}_k(s_k(t,\mathbf{r}))\cdot (\mathbf{r}-{\mathbf{p}}_k(s_k(t,\mathbf{r})))}{{\left|\mathbf{r}-{\mathbf{p}}_k(s_k(t,\mathbf{r}))\right|}^2\cdot (h_k(s_k(t,\mathbf{r}),\mathbf{r}))}\right)$$

$$\left\{\begin{array}[]{rll}\displaystyle\frac{1}{c^{2}}\partial_{t}^{2}v-\Delta v&=0,&\qquad\mbox{ in }(0,T)\times\Omega,\\ v(T,\boldsymbol{r})&=0,&\qquad\mbox{ in }\Omega,\\ \partial_{t}v(T,\boldsymbol{r})&=0,&\qquad\mbox{ in }\Omega.\end{array}\right.$$

$${\mathcal{R}}_{\varphi }(v)\equiv -{\int }_0^T{\int }_{\Gamma }\varphi (t,\mathbf{r})v(t,\mathbf{r})dS(\mathbf{r})dt,v\in \mathcal{W}.$$

$${\mathcal{R}}_{\varphi }(v)=$$

$$-{\int }_0^T{\int }_{\Omega }u(t,\mathbf{r})\left(\frac{1}{c^2}{\partial }_t^{2}v(t,\mathbf{r})-\Delta v(t,\mathbf{r})\right)dV(\mathbf{r})dt$$

$$=$$

$$f_{n,\epsilon }(t,\mathbf{r};\tau )=$$

$$g_{n,\epsilon }(t,\mathbf{r};\tau )=$$

$$h_{n,\epsilon }(t,\mathbf{r};\tau )=$$

$$n=1,2,\cdots ,$$

$$i_{n,\epsilon }(t,\mathbf{r};\tau )=$$

$$j_{n,\epsilon }(t,\mathbf{r};\tau )=$$

Full text

{amscode}

35L05, 35R30,65M32

1 Introduction

Reconstruction problem for wave sources often arises in science, engineering and medical fields, and has many important applications in these areas, e.g. reconstruction of seismic source, identification of noisy sound source in the human body, and so on [1, 2, 3]. Such kind of problem can be formulated as an inverse source problem for wave equation. In this paper, we consider a simplified case such that the media is homogeneous and isotropic, and measurement of the wave field is given on the whole boundary of the domain. Then we can formulate the problem as follows:

Let $u$ be a wave field described as a solution of the following initial-boundary value problem for three dimensional scalar wave equation:

[TABLE]

where $\boldsymbol{r}=(x,y,z)$, $\Omega\subset{\mathbb{R}}^{3}$ is a simply connected bounded domain with $C^{\infty}$-boundary $\Gamma=\partial\Omega$, $c>0$ is the given wave propagation speed, $T>2\cdot{\mathrm{diam}}\,\Omega/c>0$ is a given constant, and $F(t,\boldsymbol{r})$ describes the wave source. In (1), suppose that the wave source $F$ is unknown, and consider the problem to reconstruct unknown $F$ from observation $\phi$ on $\Gamma$ given by

[TABLE]

where $\partial_{\nu}u$ denotes the outward normal derivative of $u$ on $\Gamma$.

Many researchers discussed the reconstruction of wave sources in theoretical and numerical points of view, e.g. [4, 5, 6, 7, 8]. Unfortunately, the solution of this inverse problem is not unique in general, and hence we usually restrict the source term to some ideal models. In this paper, we assume that the source term is described by multiple moving point sources

[TABLE]

or moving dipole sources

[TABLE]

In (2), $K$ denotes the number of point sources, and $\boldsymbol{p}_{k}(t)\equiv(p_{k,x}(t),\,p_{k,y}(t),\,p_{k,z}(t))\in D$ and $q_{k}(t)\in{\mathbb{R}}$ denote the location and magnitude of $k$-th point source at $t$, respectively, where $D$ is a compact subset of $\Omega$. Note that we need not specify $D$ explicitly. The symbol $\delta$ describes the Dirac’s delta distribution which is understood as a linear functional

[TABLE]

for any $v\in C(\overline{\Omega})$. Also in (3), $K$ denotes the number of dipole sources, and $\boldsymbol{p}_{k}(t)\in D$ and $\boldsymbol{m}_{k}(t)\equiv(m_{k,x}(t),\,m_{k,y}(t),\,m_{k,z}(t))\in{\mathbb{R}}^{3}$ denote the location and dipole moment of $k$-th dipole source at $t$, respectively. The gradient of $\delta$ is understood as

[TABLE]

for any $v\in C^{1}(\Omega)\cap C(\overline{\Omega})$. Hence, we consider the solution $u$ of (1) in a weak sense, i.e. $u\in C^{1}((0,T];H^{1}_{0}(\Omega))$ that satisfies

[TABLE]

for any $v\in C^{\infty}([0,T]\times\overline{\Omega})$. The right hand side term ${\mathcal{F}}(v)$ is expressed by

[TABLE]

for moving point sources (2), and

[TABLE]

for moving dipole sources (3).

Reconstruction methods of unknown sources can be categorised into two types, one is based on an optimisation technique like the least squares method[5, 9, 10, 11, 12], and the other is based on an algebraic idea [13, 14, 15, 16]. The former method has high versatility, but it usually requires large computational cost because we need to solve partial differential equation iteratively. In the latter method, the problem is usually reduced into algebraic equations of parameters of unknown sources. This kind of method has a merit on the computational cost since we need not to solve partial differential equation iteratively instead of a drawback on relatively low versatility. In this paper, we focus our interest on the latter type of methods.

The first work on algebraic reconstruction methods for wave sources was done by El Badia and Ha Duong in 2001[13, 14]. They assumed that several point sources are fixed in $D$, i.e. $\boldsymbol{p}_{k}(t)\equiv\boldsymbol{p}_{k}\in D$, and derived a reconstruction procedure from observations $\phi$ on the whole of $\Gamma$ and the final state $u(T,\boldsymbol{r})$ and $\partial_{t}u(T,\boldsymbol{r})$ in $\Omega$. The key ideas of their method are the reciprocity gap functional [17] with respect to the space variable $\boldsymbol{r}$ and the Fourier transform with respect to the time variable $t$. Since their method needs to compute the inverse Fourier transform in the reconstruction procedure, we can not reconstruct unknown source is simultaneously on observing data. In 2011, Inui, Ohnaka and the author extended the idea of [14], and propose a new algebraic procedure in which one can reconstruct the locations and magnitudes of several point sources almost simultaneously on observations [16]. We can also find some researches which apply the reciprocity gap functional to localise a point wave source[18, 19]. For the reconstruction of moving wave source, only a few papers are published in present(e.g. [20, 21, 22, 23]), and in these papers unknown source have been assumed as a single point source. The purpose of this paper is to extend our result [16] to the case of multiple moving point and dipole sources, and give a real-time reconstruction procedure, where the phrase ’real-time’ means simultaneously on observing data.

The contents of this paper are as follows. Section 2 shows some regularity results for the solution of (1), especially, for observation data $\phi$. In section 3, we describe the detail of our reconstruction method; First, we define the reciprocity gap functional for the wave equation. Next, we give an appropriate set of functions, and derive algebraic relations between the reciprocity gap functionals for these functions and parameters of unknown sources. Based on these relations, we propose a real-time reconstruction procedure for unknown sources. In section 4, we perform some numerical experiments and discuss the effectiveness of our reconstruction procedure.

Throughout this paper, we use the following notations for some functional spaces. The Lebesgue space of the square integrable functions on $\Omega$ is denoted by $L^{2}(\Omega)$, and the set of all functions $u$ of which the weak derivative $D^{\alpha}u\in L^{2}(\Omega)$ for every multiindex $|\alpha|\leq m$ by $H^{m}(\Omega)$. $H^{m}_{0}(\Omega)$ denotes the closure of ${\mathcal{D}}(\Omega)$ in $H^{m}(\Omega)$ where ${\mathcal{D}}(\Omega)$ is the set of all infinitely differentiable functions with the compact support in $\Omega$. We denote the dual space of $H^{m}_{0}(\Omega)$ by $H^{-m}(\Omega)$.

2 Some regularity results for observation data

Before discussion of the reconstruction of wave sources, we show a regularity result on the solution of initial-boundary value problem (1) for moving point and dipole sources, especially, the regularity of the observation data $\phi=\partial_{\nu}u$ on $\Gamma$.

Proposition 2.1**.**

Let $\ell$ be a non-negative integer. For moving point sources (2), assume that $\boldsymbol{p}_{k}\in C^{\ell+2}([0,T];D),\ q_{k}\in C^{\ell+1}([0,T];\mathbb{R})$, $|d_{t}{\boldsymbol{p}}_{k}(t)|<c$ and $q_{k}(0)=d_{t}{q}_{k}(0)=d^{2}_{t}{q}_{k}(0)=\cdots=d_{t}^{(\ell+1)}q_{k}(0)=0$, where $d_{t}$ denotes the derivative with respect to $t$. Then, the solution $u$ of (1) satisfies $u\in C([0,T];L^{2}(\Omega))$ and $\partial_{t}u\in C([0,T];L^{2}(\Omega))$. Specifically, the restriction $u$ on $[0,T]\times(\Omega\backslash D)$ satisfies

[TABLE]

and the normal derivative $\partial_{\nu}u$ on $\Gamma$ satisfies

[TABLE]

For moving dipole sources (3), assume that $\boldsymbol{p}_{k}\in C^{\ell+2}([0,T];D),\ \boldsymbol{m}_{k}\in C^{\ell+2}([0,T];{\mathbb{R}}^{3})$, $|d_{t}{\boldsymbol{p}}_{k}(t)|<c$, and $\boldsymbol{m}_{k}(0)=d_{t}{\boldsymbol{m}}_{k}(0)=d_{t}^{2}{\boldsymbol{m}}_{k}(0)=\cdots=d_{t}^{(\ell+2)}\boldsymbol{m}_{k}(0)=\boldsymbol{0}$. Then, the solution $u$ satisfies $u\in C([0,T];H^{-1}(\Omega))$, and $\partial_{t}u\in C([0,T];H^{-1}(\Omega))$. Specifically, the restriction $u$ on $[0,T]\times(\Omega\backslash D)$ and the normal derivative $\partial_{\nu}u$ have the same regularity properties (7)-(10) as for moving point sources.

Proof.

First, we consider the case where the source term $F$ is expressed by moving point sources (2). Let $u_{\mathrm{N}}$ be the solution of the following initial value problem for the wave equation in the whole space ${\mathbb{R}}^{3}$:

[TABLE]

Then, $u_{\mathrm{N}}$ is explicitly expressed by the following extended Liénard-Wiechert retarded potential[24]:

[TABLE]

where $s_{k}(t,\boldsymbol{r})$ is determined as a solution $s$ of the equation

[TABLE]

for each $t,\boldsymbol{r}$ and $k$, and

[TABLE]

Since $\boldsymbol{p}_{k}\in C^{\ell+2}([0,T];D)$ and $|d_{t}{\boldsymbol{p}}_{k}(t)|<c$, $s_{k}(t,\boldsymbol{r})$ is uniquely determined and $s_{k}\in C^{\ell+1}(([0,T]\times\Omega)\backslash\bigcup_{t\in(0,T)}\{(t,\boldsymbol{p}_{k}(t))\})\cap C([0,T]\times\Omega)$, and hence we have $h_{k}\in C^{\ell+1}(([0,T]\times\Omega)\backslash\bigcup_{t\in(0,T)}\{(t,\boldsymbol{p}_{k}(t))\})$. For the derivation of (11), see e.g. Appendix A in [21]. The restriction of $u_{\mathrm{N}}$ on $[0,T]\times\Omega$ satisfies $u_{\mathrm{N}}|_{[0,T]\times\Omega}\in C([0,T];L^{2}(\Omega))$ since $s_{k}(t,\boldsymbol{p}_{k}(t))=t$ and $|\boldsymbol{r}-\boldsymbol{p}_{k}(s_{k}(t,\boldsymbol{r}))|=O(|\boldsymbol{r}-\boldsymbol{p}_{k}(t)|)$ as $\boldsymbol{r}\rightarrow\boldsymbol{p}_{k}(t)$. Moreover, the restriction to $[0,T]\times({\mathbb{R}^{3}}\backslash D)$ satisfies $u_{\mathrm{N}}|_{[0,T]\times({\mathbb{R}^{3}}\backslash D)}\in C^{\ell+1}([0,T]\times({\mathbb{R}^{3}}\backslash D))$ and $\nabla u_{\mathrm{N}}|_{[0,T]\times({\mathbb{R}^{3}}\backslash D)}\in C^{\ell}([0,T]\times({\mathbb{R}^{3}}\backslash D))$. Hence, $u_{\mathrm{N}}|_{[0,T]\times\Gamma}\in C^{\ell+1}([0,T]\times\Gamma)\subset H^{\ell+1}([0,T]\times\Gamma)$ and $\bm\nu\cdot\nabla u_{\mathrm{N}}|_{[0,T]\times\Gamma}\in C^{\ell}([0,T]\times\Gamma)\subset H^{\ell}([0,T]\times\Gamma)$ since $\Gamma$ is in $C^{\infty}$-class, where $\bm\nu$ is the unit vector normal to $\Gamma$.

Next, let $u_{\mathrm{H}}$ be the solution of initial-boundary value problem of the homogeneous wave equation with inhomogeneous Dirichlet condition on $\Gamma$:

[TABLE]

Owing to Remark 2.10 in [25], the solution $u_{\mathrm{H}}$ of (12) exists in $C([0,T];H^{\ell+1}(\Omega))$ since $u_{\mathrm{N}}|_{[0,T]\times\Gamma}$ satisfies all compatibility condition $\partial_{t}^{k}u_{\mathrm{N}}(0,\boldsymbol{r})=0$ for $k=0,1,2,\cdots,\ell+1$ on $\Gamma$. Moreover, $\partial^{m}_{t}u_{\mathrm{H}}\in C([0,T];H^{\ell-m+1}(\Omega))$ and $\partial_{\nu}u_{\mathrm{H}}\in H^{\ell}([0,T]\times\Gamma)$.

Let $u=u_{\mathrm{H}}+u_{\mathrm{N}}$, then $u$ becomes the solution of initial-boundary value problem (1) in $C([0,T];L^{2}(\Omega))$, and satisfies (7)-(10).

Next, we consider the case where the source term is expressed by dipole sources model (3). As same as for point sources model, let us consider the solution $u_{\mathrm{N}}$ of the following initial value problem in $\mathbb{R}^{3}$:

[TABLE]

Then, $u_{\mathrm{N}}$ is explicitly expressed as

[TABLE]

We can derive (13) using a similar discussion as the one for moving point sources, but we omit it. Under the assumption for $\boldsymbol{p}_{k}$ and $\boldsymbol{m}_{k}$, one can see that $u_{\mathrm{N}}\in C^{\ell+1}([0,T];H^{-1}(\Omega))$. Also for restrictions $u_{\mathrm{N}}|_{[0,T]\times({\mathbb{R}^{3}\backslash D})}$, $u_{\mathrm{N}}|_{[0,T]\times\Gamma}$ and $\bm\nu\cdot\nabla u_{\mathrm{d}}|_{[0,T]\times\Gamma}$, we can show the same regularities as for moving point sources. Then, using the same discussion, we can derive that the solution $u$ of (1) is in $C([0,T];H^{-1}(\Omega))$, and (7)-(10).

∎

3 Reconstruction method

3.1 Reciprocity gap functional

In this section, we present a reconstruction method for moving point and dipole sources from observation $\phi=\partial_{\nu}u$ on $\Gamma$. The key idea of our method is the reciprocity gap functional that is defined on the subspace of test functions $v$ in (4). This idea is widely applied to various inverse problems, e.g. inverse conductivity problems and inverse scattering problems [26, 17, 27, 28, 29, 30, 15, 31]. First, we show a definition of the reciprocity gap functional for scalar wave equations.

Let ${\mathcal{W}}\subset C^{\infty}([0,T]\times\overline{\Omega};{\mathbb{C}})$ be a set of complex-valued functions $v$ that satisfy the homogeneous wave equation with vanishing condition at $t=T$:

[TABLE]

For given observation data $\phi\in L^{2}((0,T)\times\Gamma)$, we define the reciprocity gap functional $\mathcal{R}_{\phi}$ on ${\mathcal{W}}$ as follows:

[TABLE]

Then, since $\phi=\partial_{\nu}u$ and $u$ satisfies the weak form (4), we obtain

[TABLE]

The equation (15) gives the relation between the reciprocity gap functional ${\mathcal{R}}_{\phi}(v)$ and the source term $F$, and suggests that we can reconstruct the unknown source term from ${\mathcal{R}}_{\phi}(v)$ for appropriate choice of functions $v$. In the following subsections, we give our choice of functions $v\in\cal W$, and propose reconstruction procedures of moving point sources and dipole sources.

3.2 Reconstruction of moving point sources

Throughout this subsection, we assume that $\boldsymbol{p}_{k}\in C^{6}([0,T];D),\ q_{k}\in C^{5}([0,T];\mathbb{R})$ and $q_{k}(0)=d_{t}{q}_{k}(0)=d_{t}^{2}{q}_{k}(0)=\cdots=d_{t}^{5}q_{k}(0)=0$. Then, from Proposition 2.1 and since $[0,T]\times\Gamma$ is a three dimensional smooth manifold, observation data $\phi$ is in $H^{4}((0,T)\times\Gamma)\subset C^{2}((0,T)\times\Gamma)$. For the reconstruction of moving point sources, we choose the following five sequences of functions in $\mathcal{W}$:

[TABLE]

where $\tau\in I_{\varepsilon}\equiv\left[\max\left\{0,\sup_{\boldsymbol{r}\in\Omega}z/c+\varepsilon\right\},\,\min\left\{T,T+\inf_{\boldsymbol{r}\in\Omega}z/c-\varepsilon\right\}\right]\subset[0,T]$, $0<\varepsilon\ll 1$, and $\eta_{\varepsilon}\in C^{\infty}({\mathbb{R}};{\mathbb{R}})$ denotes the standard mollifier function with the support $[-\varepsilon,\varepsilon]$ (e.g. Appendix C in [32]). We note that sequences $\{f_{n,\varepsilon}\},\ \{g_{n,\varepsilon}\}$ and $\{h_{n,\varepsilon}\}$ have already applied to the reconstruction of fixed point sources [16]. Supplemental sequences $\{i_{n,\varepsilon}\}$ and $\{j_{n,\varepsilon}\}$ are used to treat the effect of moving velocities of sources.

Due to the assumptions for $\boldsymbol{p}_{k}$ and $q_{k}$, the observation data $\phi$ is in $C^{2}([0,T]\times\Gamma)$. Then reciprocity gap functionals ${\mathcal{R}}_{\phi}(f_{n,\varepsilon}),\ {\mathcal{R}}_{\phi}(g_{n,\varepsilon}),\cdots,{\mathcal{R}}_{\phi}(j_{n,\varepsilon})$ converge as $\varepsilon\rightarrow+0$. Now, let

[TABLE]

Then we establish

[TABLE]

and similarly,

[TABLE]

since the observation data $\phi$ is in $C^{2}((0,T)\times\Gamma)$, where $d_{\tau}$ denotes the derivative with respect to $\tau$.

Also since ${\mathcal{F}}(f_{{n,\varepsilon}}),\ {\mathcal{F}}(g_{{n,\varepsilon}}),\cdots,{\mathcal{F}}(i_{{n,\varepsilon}})$ converge as $\varepsilon\rightarrow+0$, we obtain the following formulae that express explicit relations to the parameters of moving point sources:

[TABLE]

[TABLE]

where

• •

$p_{k,xy}(t_{k}(\tau))=p_{k,x}(t_{k}(\tau))+\mathrm{i}p_{k,y}(t_{k}(\tau))$,

• •

$t_{k}(\tau)$ is the unique solution $t$ of the equation

[TABLE]

for each $k$ and $\tau$,

• •

$\xi_{k}(\tau)$ is the derivative of $t_{k}(\tau)$ and can be expressed by

[TABLE]

• •

$\hat{R}_{h_{n}},\ \hat{R}_{i_{n}}$ and $\hat{R}_{j_{n}}$ are defined by

[TABLE]

[TABLE]

In (27)-(30), we omit the argument $(t_{k}(\tau))$ on $q_{k},\ p_{k,xy},\ p_{k,z}$ and their derivatives, and the argument $(\tau)$ on $\xi_{k}$ to simplify the expression, e.g. $q_{k}(t_{k}(\tau))$ is simplified as $q_{k}$. In the following, we use these notations if we do not need to show these arguments explicitly. Derivations of (26)-(30) are given in Appendix A. We note that equations (26)-(28) are similar to the results for fixed point sources [16], but some supplemental terms arise due to the effect of moving velocities of sources, e.g. $\xi_{k}(\tau)$ and ${\mathrm{d}}_{\tau}(p_{k,xy}(t_{k}(\tau)))$.

Note 3.1. Since functions $f_{n,\varepsilon},\ g_{n,\varepsilon},\cdots,\ j_{n,\varepsilon}$ have propagation property along only $z$-axis, the reciprocity gap functionals of these functions have ‘retarded’ property only on $z$-coordinate. Also we can obtain an another expression of the perturbation factor $\xi_{k}(\tau)$ due to the $z$-component of the velocity of point sources:

[TABLE]

Equation (36) shows that the perturbation factor $\xi_{k}$ effects similarly as the factor $1/h$ in the extended Liénard-Wiechert retarded potential (11) and (13).

Using expressions (26)-(30), we obtain the following reconstruction theorem for moving point sources:

Theorem 3.1**.**

Let $I_{0}\equiv\left[\max\{0,\sup_{\boldsymbol{r}\in\Omega}z/c\},\min\{T+\inf_{\boldsymbol{r}\in\Omega}z/c,T\}\right]\subset[0,T]$. For each $\tau\in I_{0}$, let $K(\tau)$ be the number of point sources such that $q_{k}(t_{k}(\tau))\neq 0$. Assume that $K(\tau)\leq K_{M}$ for given $K_{M}$, and $p_{j,xy}(t_{j}(\tau))\neq p_{k,xy}(t_{k}(\tau))$ if $j\neq k$. Then, we can identify $K(\tau)$ from the reciprocity gap functionals ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=0,1,2,\cdots,2K_{M}$. Also we can uniquely reconstruct $\boldsymbol{p}_{k}(t_{k}(\tau))$ and $q_{k}(t_{k}(\tau))$, $k=1,2,\cdots,K(\tau)$ from

• •

${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$,

• •

${\mathcal{R}}_{\phi}(g_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$,

• •

${\mathcal{R}}_{\phi}(h_{n})(\tau),\ n=1,2,3,\cdots,K(\tau)$,

• •

${\mathcal{R}}_{\phi}(i_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$,

• •

${\mathcal{R}}_{\phi}(j_{n})(\tau),\ n=1,2,3,\cdots,K(\tau)$.

Proof.

We show the proof of the theorem in the following five steps which describe the procedure of the reconstruction process.

Step 1.

Identify $K(\tau)$ from ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=0,1,2,\cdots,2K_{M}$.

Step 2.

Reconstruct $p_{k,xy}(t_{k}(\tau)),\ k=1,2,\cdots,K(\tau)$, and identify perturbed magnitudes $q_{k}(t_{k}(\tau))\xi_{k}(\tau)$ from ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$.

Step 3.

Reconstruct $p_{k,z}(t_{k}(\tau)),\ k=1,2,\cdots,K(\tau)$ from ${\mathcal{R}}_{\phi}(g_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$ and ${\mathcal{R}}_{\phi}(h_{n})(\tau),\ n=1,2,3,\cdots,K(\tau)$.

Step 4.

Identify $d_{\tau}(p_{k,z}(t_{k}(\tau))),\ k=1,2,\cdots,K(\tau)$ from ${\mathcal{R}}_{\phi}(i_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$ and ${\mathcal{R}}_{\phi}(j_{n})(\tau),\ n=1,2,3,\cdots,K(\tau)$.

Step 5.

Compute $\xi_{k}(\tau)$ using $d_{\tau}({p}_{k,z}(t_{k}(\tau)))$ for each $k$, and reconstruct the magnitude $q_{k}(t_{k}(\tau))$ from its perturbed value $q_{k}(t_{k}(\tau))\xi_{k}(\tau)$.

We show the detail of each step bellow.

Step 1. Define $L\times L$ Hankel matrix

[TABLE]

Then, from (26) and using corollary 3 in [31], we can determine $K(\tau)$ by

[TABLE]

Equation (37) shows that $\det H_{L,0}(\tau)$ vanishes for $L\geq K(\tau)+1$. Hence we can determine $K(\tau)$ from ${\mathcal{R}}_{\phi}(f_{n})(\tau)$ for $n=0,1,2,\cdots,2K_{M}$ if $K(\tau)\leq K_{M}$

Step 2. From the definition of $H_{L,\mu}(\tau)$ and using Theorem 2 in [29], we can reconstruct $p_{k,xy}(t_{k}(\tau)),\ k=1,2,\cdots,K(\tau)$ as eigenvalues of the matrix $(H_{K(\tau),0}(\tau))^{-1}H_{K(\tau),1}(\tau)$. For this computation, we need $\mathcal{R}_{\phi}(f_{n})(\tau)$ for $n=0,1,2,\cdots,2K(\tau)-1$.

Now, let us define $K(\tau)\times K(\tau)$-matrix

[TABLE]

and ${\bm\zeta}(\tau),{\bm\beta}(\tau)\in{\mathbb{C}}^{K(\tau)}$ by

[TABLE]

Then, equations (26) for $n=0,1,2,\cdots,K(\tau)-1$ are rewritten by

[TABLE]

Since $V(\tau)$ is a Vandermonde type matrix and from the assumption that $p_{j,xy}(t_{j}(\tau))\neq p_{k,xy}(t_{k}(\tau))$ for $j\neq k$, we derive $\det V(\tau)\neq 0$ . Then, the equation (38) is uniquely solvable, and we can identify the perturbed magnitudes $q_{k}(t_{k}(\tau))\,\xi_{k}(\tau),\ k=1,2,\cdots,K(\tau)$.

Step 3. Owing to equation (28), if we know values of ${\mathcal{R}}_{\phi}(h_{n})(\tau)$ and $\hat{R}_{h_{n}}$ for $n=1,2,\cdots,K(\tau)$, we can identify $q_{k}\,\xi_{k}\cdot p_{k,z},\ k=1,2,\cdots,K(\tau)$ as the unique solution of the system of linear equation:

[TABLE]

Then, dividing each solution $q_{k}\xi_{k}\cdot p_{k,z}$ of (39) by the perturbed magnitude $q_{k}\xi_{k}$ that is identified in step 2, we obtain $p_{k,z}$. Hence we consider the estimation of $\hat{R}_{h_{n}}$.

In the expression (33) of $\hat{R}_{h_{n}}$, $d_{\tau}(q_{k}\xi_{k})$ and $d_{\tau}(p_{k,xy})$ remain unknowns. We apply (27) to identify these unknowns. Let us define ${\bm\psi}_{k},\ {\bm\psi}^{\prime}_{k}\in{\mathbb{C}}^{2K(\tau)}$ by

[TABLE]

and the matrix $\tilde{V}(\tau)\in{\mathbb{C}}^{2K(\tau)\times 2K(\tau)}$ by

[TABLE]

Also define ${\bm\iota}(\tau),\ {\bm\gamma}(\tau)\in{\mathbb{C}}^{2K(\tau)}$ by

[TABLE]

Then, equations (27) for $n=0,1,2,\cdots,2K(\tau)-1$ are rewritten as

[TABLE]

The matrix $\tilde{V}(\tau)$ is a generalised Vandermonde-type matrix, and we establish the following equality:

[TABLE]

A proof of (42) is given in Appendix B. From the assumption that $p_{j,xy}(t_{j}(\tau))\neq p_{k,xy}(t_{k}(\tau))$ if $j\neq k$, we obtain $\det\tilde{V}(\tau)\neq 0$. Then, the equation (41) is uniquely solvable, and we can identify $d_{\tau}(q_{k}\,\xi_{k})$ and $d_{\tau}(p_{k,xy}),\ k=1,2,\cdots,K(\tau)$. This concludes that we can reconstruct $p_{k,z}(t_{k}(\tau)),\ k=1,2,\cdots,K(\tau)$.

Step 4. Owing to (30), we can identify $q_{k}\xi_{k}\cdot d_{\tau}(p_{k,z})$ as the unique solution of the following system of linear equation if we know $\mathcal{R}_{\phi}(j_{n})(\tau)$ and $\hat{R}_{j_{n}}$ for $n=1,2,\cdots,K(\tau)$:

[TABLE]

From the expression (35) of $\hat{R}_{j_{n}}$, we can see that $d_{\tau}^{2}(q_{k}\xi_{k})$ and $d^{2}_{\tau}(p_{k,xy})$ are unknown in the right hand side term of (43). Using the same idea as Step 3, these unknowns can be identified from ${\mathcal{R}}_{\phi}(i_{n})(\tau),\ n=0,1,2,\cdots,2K(\tau)-1$ since ${\mathcal{R}}_{\phi}(i_{n})(\tau)$ and $\hat{R}_{i_{n}}$ are given by (29) and (34), respectively. Then, we have identified all unknowns in the right hand side terms of (43), and we can reconstruct $d_{\tau}(p_{k,z}),\ k=1,2,\cdots,K(\tau)$ similarly as Step 3.

Step 5. Since we have identified $d_{\tau}(p_{k,z})$ in Step 4, we can compute the perturbation term $\xi_{k}$ by (32), and reconstruct magnitude $q_{k}$ of each point source from perturbed magnitude $q_{k}\,\xi_{k}$ for $k=1,2,\cdots,K(\tau)$.

Then, we have reconstructed all parameters of moving point sources at $t=t_{k}(\tau),\ k=1,2,\cdots,K(\tau)$, and this completes the proof of Theorem 3.1.

∎

Note 3.2. We can use alternative functions that propagate to positive direction with respect to $z$-coordinate, i.e.

[TABLE]

In this case, $t_{k}(\tau)$ is replaced by the solution of $t$ of the equation

[TABLE]

and $\xi_{k}(\tau)$ is replaced by

[TABLE]

Note 3.3. In practical cases, it is difficult to apply the condition (37) to the identification of $K(\tau)$ because of noise in the observation data and errors of numerical computation. In section 4, we propose another heuristic criterion to identify the number of sources.

Note 3.4. In our reconstruction procedure, we reconstruct the locations and magnitude of points sources at time $t_{k}(\tau)$ for each $k$ and $\tau$. We can easily estimate $t_{k}(\tau)$ for each $k$ and $\tau$ by replacing $t$ of (31) by $t_{k}(\tau)$, i.e.

[TABLE]

Here, we do not need to solve this equation with respect to $t_{k}(\tau)$ but only substitute $p_{k,z}(t_{k}(\tau))$ that was identified in Step 3.

Note 3.5. Since we have obtained $p_{k,z}(t_{k}(\tau))$ in Step 3, we may compute $d_{\tau}p_{k,z}(t_{k}(\tau))$ using numerical differentiation in stead of Step 4. (Here, we call Step 4’.) We compare these two methods in section 4.

Note 3.6. The merits of our reconstruction procedure are (a) real-timeness and (b) independence on the initial condition. Reconstruction methods based on an optimisation procedure such as the least squares method usually need whole time observation data and the initial condition because we need to solve the initial-boundary value problem. On the other hand, the proof of Theorem 3.1 shows that our method needs observations in the interval $[\tau-\sup_{\boldsymbol{r}\in\Omega}z/c,\ \tau-\inf_{\boldsymbol{r}\in\Omega}z/c]$ to reconstruct unknown parameters at $t=t_{k}(\tau)$. This shows the real-timeness and independence on the initial condition of our reconstruction procedure.

3.3 Reconstruction of moving dipole sources

Next, we consider the reconstruction of moving dipole sources. Here, we assume that $\boldsymbol{p}_{k}\in C^{6}([0,T];D),\ \boldsymbol{m}_{k}\in C^{6}([0,T];\mathbb{R}^{3})$ and $\boldsymbol{m}_{k}(0)=d_{t}{\boldsymbol{m}}_{k}(0)=d_{t}^{2}{\boldsymbol{m}}_{k}(0)=\cdots=d_{t}^{6}\boldsymbol{m}_{k}(0)=\boldsymbol{0}$, then the observation data $\phi$ is in $H^{4}([0,T]\times\Gamma)\subset C^{2}([0,T]\times\Gamma)$. We also add an assumption that $m_{k,z}(t)\equiv 0$ for all $k$, i.e. the moment of each dipole source is expressed by $\boldsymbol{m}_{k}(t)=(m_{k,x}(t),\,m_{k,y}(t),0)$. For the reconstruction of moving dipole sources, we apply the same five sequences of functions $\{f_{n,\varepsilon}\},\{g_{n,\varepsilon}\},\cdots,\{j_{n,\varepsilon}\}$ as for the reconstruction of moving point sources. Then, we obtain the following expressions of ${\mathcal{F}}(f_{n})(\tau)$, ${\mathcal{F}}(g_{n})(\tau)$, $\cdots$, ${\mathcal{F}}(j_{n})(\tau)$ using the parameters of moving dipole sources:

[TABLE]

where

• •

$m_{k,xy}\equiv m_{k,x}(t_{k}(\tau))+\mathrm{i}m_{k,y}(t_{k}(\tau))$,

• •

$\hat{R}_{h_{n}}$, $\hat{R}_{i_{n}}$, and $\hat{R}_{j_{n}}$ are expressed by

[TABLE]

[TABLE]

• •

$t_{k}(\tau)$ and $\xi_{k}(\tau)$ are defined by (31) and (32), respectively.

Derivations of equations (44)-(48) are also given in Appendix A.

Similarly to the case for moving point sources, we can establish the following reconstruction theorem for moving dipole sources:

Theorem 3.2**.**

Let $I_{0}$ be the same interval as in Theorem 3.1. For each $\tau\in I_{0}$, let $K(\tau)$ be the number of dipole sources such that $\boldsymbol{m}_{k}(t_{k}(\tau))\neq\boldsymbol{0}$. Assume that $K(\tau)\leq K_{M}$ for given $K_{M}$, $p_{j,xy}(t_{j}(\tau))\neq p_{k,xy}(t_{k}(\tau))$ if $j\neq k$, and $m_{k,z}(t)\equiv 0$. Then, we can identify $K(\tau)$ from the reciprocity gap functionals ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=1,2,\cdots,2K_{M}+1$. Also we can uniquely reconstruct $\boldsymbol{p}_{k}(t_{k}(\tau))$ and $\boldsymbol{m}_{k}(t_{k}(\tau))$ for $k=1,2,\cdots,K(\tau)$ from

• •

${\mathcal{R}}(f_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$,

• •

${\mathcal{R}}(g_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$,

• •

${\mathcal{R}}(h_{n})(\tau),\ n=2,3,4,\cdots,K(\tau)+1$,

• •

${\mathcal{R}}(i_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$,

• •

${\mathcal{R}}(j_{n})(\tau),\ n=2,3,4,\cdots,K(\tau)+1$.

Proof.

Comparing expressions (44)-(48) to (26)-(30), we can find similar relations between the reciprocity gap functionals and parameters of dipole sources as for moving point sources. Each step of the reconstruction procedure in the proof of Theorem 3.1 is modified as follows:

Steps 1 and 2. Dividing (44) by $n$, we have

[TABLE]

The right hand side term of (49) has the same form as (26) by replacing $q_{k}(t_{k}(\tau))$ by $m_{k,xy}(t_{k}(\tau))$ and the power $n$ by $n-1$. Hence, we can identify $K(\tau)$ using $H_{k,1}(\tau)$ computed from ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=1,2,\cdots,2K_{M}+1$, and reconstruct $p_{k,xy}(t_{k}(\tau))$, and $m_{k,xy}(t_{k}(\tau))\xi_{k}(\tau)$ from ${\mathcal{R}}_{\phi}(f_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$ using the same procedure as in Steps 1 and 2 of the proof of Theorem 3.1.

Steps 3 and 4. Similarly to Steps 1 and 2, by dividing (45)-(48) by $n$, we can obtain the same forms as (27)-(30) by replacing $q_{k}(t_{k}(\tau))$ by $m_{k,xy}(t_{k}(\tau))$ and the power $n$ by $n-1$. Hence, we can identify $p_{k,z}(t_{k}(\tau))$ from ${\mathcal{R}}_{\phi}(g_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$ and ${\mathcal{R}}_{\phi}(h_{n})(\tau),\ n=2,\cdots,K(\tau)+1$, and $d_{\tau}(p_{k,z}(t_{k}(\tau)))$ from ${\mathcal{R}}_{\phi}(i_{n})(\tau),\ n=1,2,3,\cdots,2K(\tau)$ and ${\mathcal{R}}_{\phi}(j_{n})(\tau),\ n=2,\cdots,K(\tau)+1$ using the same procedure as in Steps 3 and 4 of the proof of Theorem 3.1.

Step 5. Finally, we can compute $\xi_{k}(\tau)$ with (32) using $d_{\tau}(p_{k,z}(t_{k}(\tau)))$ identified in Step 4, and reconstruct each moment $m_{k,xy}(t_{k}(\tau))$ from its perturbed value $m_{k,xy}(t_{k}(\tau))\,\xi_{k}(\tau)$. Then, we have reconstructed all parameters of moving dipole sources.

∎

Note 3.7. If we remove the assumption $m_{k,z}(t)\equiv 0$, the expression (44) becomes

[TABLE]

Then, the idea in Step 1 can not be applied. Reconstruction method for such cases remains an open problem.

4 Numerical Experiments

We now present some numerical experiments for the reconstruction procedure proposed in section 3. Throughout the experiments, we set $\Omega=\{\boldsymbol{r}\,|\,|\boldsymbol{r}|\leq 2\}\subset\mathbb{R}^{3}$, $c=1$ and $T=70$. To give the observation data $\phi=\partial_{\nu}u$ on $\Gamma$, we solve the initial-boundary value problem (1) numerically using the following boundary integral expression:

[TABLE]

where $G(t,\boldsymbol{r};s,\bm\rho)$ is the fundamental solution of the three-dimensional scalar wave equation defined by

[TABLE]

$\partial_{\nu(\boldsymbol{r})}$ denotes the outward normal derivative on $\Gamma$ with respect to the variable $\boldsymbol{r}$, and $u_{\mathrm{N}}$ is given by (11) for moving point sources, and (13) for moving dipole sources. The observation points of $\phi$ are arranged as

[TABLE]

where $\cos\theta_{j}$ is the $j$-th collocation point of Gauss-Legendre quadrature, and $\varphi_{k}={2\pi(k-1)}/{K}$. In the experiments, we assign $J=18,\ K=36$, and therefore 648 observation points are arranged on $\Gamma$. We give observation data at $t=t_{\ell}=\Delta t\times\ell,\ \ell=0,1,2,\cdots,T/\Delta t$ on $[0,T]$ where $\Delta t=0.1$. To simulate practical observation conditions, we add Gaussian noise to the numerical solution at each time $t=t_{\ell}$ such that

[TABLE]

where $\phi^{\mathrm{obs}}$ denotes the observation data perturbed by Gaussian noise. The application of our reconstruction procedure is performed for every $\tau=\tau_{\ell}=\Delta\tau\times\ell,\ \ell=0,1,2,\cdots,T_{e}/\Delta\tau$, where $\Delta\tau=0.1$, $T_{e}=60.0$. We approximate the surface integral on $\Gamma$ using the Gauss-Legendre quadrature with respect to $\theta$ and the trapezoidal rule with respect to $\varphi$. To approximate the derivative with respect to $\tau$, we apply the central difference, e.g. ${\mathcal{R}}_{\phi}(g_{n})(\tau)$ is approximated by

[TABLE]

4.1 Reconstruction of moving point sources

We show numerical experiments for the reconstruction of moving point sources. We arrange the parameters of three point sources as follows:

source 1.

${\boldsymbol{p}}_{1}(t)=(0.8,\ -0.3,\ 0.8\cos(0.4t)-0.2)$,

[TABLE]

source 2.

${\boldsymbol{p}}_{2}(t)=(p_{2,x}(t),\ p_{2,y}(t),\ p_{2,z}(t))$ where

[TABLE]

[TABLE]

source 3.

${\boldsymbol{p}}_{3}(t)=(p_{3,x}(t),\ p_{3,y}(t),\ p_{3,z}(t))$ where

[TABLE]

[TABLE]

Here, we set

[TABLE]

Therefore, $K(\tau)$ changes as

[TABLE]

We note that $|\boldsymbol{p}_{k}(t)|<1.6$ for any $k$ and $t$, and the range of the moving speeds of sources are $0\leq|\dot{\boldsymbol{p}}_{1}(t)|\leq 0.32$, $0.165\leq|\dot{\boldsymbol{p}}_{2}(t)|\leq 0.425$ and $0.045\leq|\dot{\boldsymbol{p}}_{3}(t)|\leq 0.333$. We also note that $q_{k}\in C^{4}([0,70]),\ k=1,2,3$.

Before we show the reconstruction results of moving point source, we compare the behaviours of magnitudes $q_{k}(t)$ with respect to $t$ and perturbed magnitudes $q_{k}(t_{k}(\tau))\xi_{k}(t_{k}(\tau))$ with respect to $\tau$ in Figure 1. We can find large perturbation on magnitudes due to the $z$-component of the moving velocities of sources, e.g. around $\tau=21,\ 30,\ 35$, and $42$ for source 1 and $30\leq\tau\leq 37$ for source 3. One of the purpose of this paper is to give a method to correct these perturbations.

Firstly, we discuss the identification of the number $K(\tau)$ of point sources. Figure 2(a) shows the behaviour of $|\det H_{k,0}(\tau)|$ in $0\leq\tau\leq T_{e}$ for $k=1,2,3,4$, and Figure 2(b) gives the behaviour of $|\det H_{k,0}(\tau)|$ at $\tau=10,\ 20$ and $40$ for observation data without noise. As we mentioned in Note 3.3, the results of Figure 2 show that we can not apply the condition (37) as it is for the identification of $K(\tau)$. However, we can observe that $|\det H_{k,0}(\tau)|$ has relatively large gap between $k=K(\tau)$ and $K(\tau)+1$. From this observation, we apply the following algorithm to identify $K(\tau)$:

Algorithm 4.1

Step 0.

Set positive parameters $\varepsilon_{0}$ and $\varepsilon_{G}$.

Step 1.

If $|\det H_{1,0}(\tau)|<\varepsilon_{0}$ and $|\det H_{1,0}(\tau)|>|\det H_{2,0}(\tau)|$, then $K(\tau)=0$.

Step 2.

If $|\det H_{1,0}(\tau)|>\varepsilon_{0}$ or $|\det H_{1,0}(\tau)|<|\det H_{2,0}(\tau)|$, then find largest number $k$ such that $|\det H_{k,0}(\tau)|/|\det H_{k-1,0}(\tau)|>\varepsilon_{G}$, and identify $K(\tau)=k$.

Figure 3 displays the identification results of $K(\tau)$ from observation data without noise and with 0.5% noise. Here, we set $\varepsilon_{0}=1.0\times 10^{-4}$ and $\varepsilon_{G}=2.5\times 10^{-2}$. From the result of Figure 3, we consider that Algorithm 4.1 works well even for observations with 0.5% noise. One can observe that identified $K(\tau)$ is smaller than the actual value around $\tau=12$ and $\tau=55$. We discuss a reason of these bad estimates later.

Now, we show the reconstruction results of unknown parameters of moving point sources. Figures 4-7 display the reconstruction results of locations and magnitudes for observation data without noise, and with 0.1%, 0.5% and 1.0% noise. In these figures, we use the time variable $t$ for the horizontal axis by plotting $(t_{k}(\tau),\boldsymbol{p}_{k}(t_{k}(\tau)))$ and $(t_{k}(\tau),q_{k}(t_{k}(\tau)))$. In Table 4.1, we show the average errors of estimated locations and magnitudes for observations without noise, and with 0.1%, 0.5%, 1.0%, and 5.0% noise. Here we define

[TABLE]

where $\widehat{\boldsymbol{p}}_{k}$ and $\widehat{q}_{k}$ denote the estimated location and magnitude of $k$-th point source, respectively, $I$ denotes the interval in which we evaluate the average error, and $|I|$ denotes the measure of $I$.

From results of Figures 4-7 and Table 4.1, we consider that our method works well if the noise of observation data is smaller than 0.5%. We can observe that if noise exceeds 1%, the influence of noise becomes unignorable in reconstruction results. Especially, estimations of magnitudes are highly affected by the observation noise. The reconstruction results become unreliable under the noise larger than 5%.

We can find that the reconstruction of location becomes bad for the source with small magnitude, e.g. the result of source 2 around $t=12$ in Figures 6(b) and 7(b), but we consider that these results are reasonable. We can also find bad reconstruction results due to small magnitude for source 2 in the interval $54.0\leq\tau<60.0$ for the observation without noise and 1.0% noise in Table 4.1, and bad identification of $K(\tau)$ around $\tau=12$ and $\tau=53$ in Figure 3(b).

Note 4.1. In Table 4.1, the error of $\boldsymbol{p}_{2}$ in the interval $54.0\leq\tau\leq 60.0$ is concentrated in $57.6\leq\tau\leq 57.8$ where $q_{2}(t_{2}(\tau))\simeq 0$. Except of these bad estimates, the average error of $\boldsymbol{p}_{2}$ in $54.0\leq\tau<60.0$ becomes $1.56E-2$ for the case without noise and $1.45E-2$ for the case with 0.1% noise.

In Figures 6 and 7, we can also find bad reconstruction results even for large magnitudes of point sources, e.g. results for sources 2 and 3 at around $t=30$. To discuss the reason of these bad estimates, we display the behaviour of distances between each two point sources in Figure 8. The results of Figures 6-8 suggest that if the arangement of point sources is well separated, then our method works well and gives good reconstruction results, however, as arrangement of point sources becomes closer, then it becomes more difficult to distinguish these sources, and our method gives bad estimations.

In Figure 3(b), we can find many wrong identification results of $K(\tau)$, especially, in $25\leq\tau\leq 38$ where $K(\tau)$ is identified as 4. However, the absolute value of magnitude of reconstructed 4th source is very small ($<0.01$) comparing to other sources, and so the 4th source can be recognised as ‘noises’ or ‘ghosts’.

Note 4.2. Since our reconstruction procedure works at every moment $\tau$, we have to identify whether the point source reconstructed at the moment $\tau$ belongs to which source reconstructed at the previous moment $\tau-\Delta\tau$, especially, when $K(\tau)$ changes. To solve this problem, we measure distances between the reconstructed sources at moments $\tau$ and $\tau-\Delta\tau$, and identify as the same source which has the minimum distance.

Throughout of the numerical results, we can observe that the observation noise affects to the estimation of magnitudes more heavily than the one of locations. This phenomena is mainly caused by the estimation of perturbation term $\xi_{k}(t_{k}(\tau))$. As an alternative method, we apply Step 4’ in Note 3.5 to obtain $dp_{k,z}(t_{k}(\tau))/d\tau$, and compare to results by Step 4. Here, we apply the central differences

[TABLE]

to estimate $dp_{k,z}(t_{k}(\tau))/d\tau$, where $p_{k,z}(t_{k}(\tau))$ is obtained in Step 3. Table 4.1 shows the average errors of estimated magnitudes using Step 4’ instead of Step 4. Comparing Tables 4.1 and 4.1, the errors are almost the same in both methods for the case where the observation noise is smaller than 0.5%. Hence we consider that both method can be used as an alternative if the noise is small.

4.2 Reconstruction of moving dipole sources

Next, we demonstrate numerical experiments for reconstruction of moving dipole sources. We consider the case where three dipole sources move in the domain $\Omega$. The locations of dipoles are arranged as same as the previous example, and moments of dipoles change as follows:

source 1.

${\boldsymbol{m}}_{1}(t)=(m_{1,x}(t),\ m_{1,y}(t),\ 0)$ where

[TABLE]

source 2.

${\boldsymbol{m}}_{2}(t)=(m_{2,x}(t),\ m_{2,y}(t),\ 0)$ where

[TABLE]

source 3.

${\boldsymbol{m}}_{3}(t)=(m_{3,x}(t),\ m_{3,y}(t),\ 0)$ where

[TABLE]

Hence $K(\tau)$ changes as

[TABLE]

We note that $\boldsymbol{m}_{k}\in C^{4}([0,70]),\ k=1,2,3$.

We first discuss the identification of the number $K(\tau)$ of dipoles. Figure 9 shows the behaviour of $|\det H_{k,1}(\tau)|$ for $0\leq\tau\leq T_{e}$. The behaviour of $|\det H_{k,1}(\tau)|$ is quite similar as for the one of $|\det H_{k,0}(\tau)|$ for point sources in Figure 2, and hence we apply Algorithm 4.1 for the identification of $K(\tau)$ using $H_{k,1}(\tau)$ instead of $H_{k,0}(\tau)$. Figure 10 displays the identification results of $K(\tau)$ from observations without noise and with 0.5% noise. Here, we set the parameters $\varepsilon_{0}=1.0\times 10^{-4}$ and $\varepsilon_{G}=2.5\times 10^{-2}$ as same as for moving point sources. From the result of Figure 10, we consider Algorithm 4.1 also works well for the identification of the number of moving dipole sources.

Next, we show the reconstruction results for locations and moments of dipole sources. Figures 11-14 display the reconstruction results from observations without noise, and with 0.1%, 0.5% and 1.0% noise. Also in Table 4.2, we show the average errors of estimated locations and moments. Here we define the average error of each moment by

[TABLE]

where $\widehat{\boldsymbol{m}}_{k}$ denotes the estimated moment of $k$-th dipole source. Same as for the reconstruction of point sources, our method works well if the noise of observation data is smaller than 0.5%. The noise of observation becomes 1%, the influence of noise becomes unignorable, and our method can not give reliable estimates under 5% noise. We can see that observation noises affect to the estimations of moments more heavily than the reconstruction of locations. The reconstruction result also becomes bad when the norm of moment is small, and when the arrangement of dipoles is clustered.

Finally, we show the average of error of estimated moments by using Step 4’ instead of Step 4 in Table 4.2. As same as for moving point source, the errors of estimation results are almost the same in both methods if observation noise is smaller than 1.0%.

5 Conclusions

We discuss a reconstruction method for multiple moving point/dipole wave sources from boundary measurements. We derive algebraic relations between the parameters of unknown sources and the reciprocity gap functionals for five sequences of functions. Based on these algebraic relations, we give a real-time reconstruction procedure for parameters of unknown sources. We examine our reconstruction procedure by numerical experiments. Numerical results show that our procedure gives reliable reconstruction results for both moving point and moving dipole sources when the observation noise is smaller than 0.5%. However, the noise becomes 1%, its influence becomes uninnorable, and the reconstruction results become unreliable under 5% noise. Such bad influences can be found more heavily on the estimation results of magnitudes/moments of sources.

We need further discussions for more complicated cases, for example, limited aperture cases, the case where the source term contains both moving point and dipole sources, and so on.

Acknowledgements

The author is grateful to Professor E. Nakaguchi in Tokyo Medical and Dental University for his useful discussions and helpful comments. Also the author would like to acknowledge fruitful the discussions at 2017 IMI Joint Use Research Program Workshop (II) ”Practical inverse problems based on interdisciplinary and industry-academia collaboration”. This work was supported in part by A3 Foresight Program “Modeling and Computation of Applied Inverse Problems” and Grant-in-Aid for Scientific Research (S) 15H05740 of Japan Society for the Promotion of Science.

Appendix A. Derivation of

(26)-(30) and (44)-(48)

Firstly, we show the derivation of (26)-(30). From the definition of $f_{n,\varepsilon}$, we have

[TABLE]

Let $\displaystyle s=t+\frac{p_{k,z}(t)}{c}$ for each $k$, then

[TABLE]

where $t_{k}(s)$ is defined in section 3.2, and we obtain

[TABLE]

where $s_{0,k}=p_{k,z}(0)/c$ and $s_{1,k}=T+p_{k,z}(T)/c$. Since $\eta_{\varepsilon}$ is the standard mollifier function, we establish (26).

Next, from the definition of $g_{n,\varepsilon}$, we have

[TABLE]

Using integration by parts, it follows that

[TABLE]

Here, we use $t^{\prime}_{k}(s)=\xi_{k}(s),\ t_{k}(s_{0,k})=0,\ t_{k}(s_{1,k})=T$, and $\eta_{\varepsilon}(s_{0,k}-\tau)=\eta_{\varepsilon}(s_{1,k}-\tau)=0$. Then, using similar computation as (53), we establish (27).

Also for ${\mathcal{F}}(h_{n})(\tau),\ {\mathcal{F}}(i_{n})(\tau)$ and ${\mathcal{F}}(j_{n})(\tau)$, we can obtain the following form by changing variable $t$ to $s$ and integration by parts:

[TABLE]

[TABLE]

Then, using similar computation as (54), we establish (28)-(30).

Next, we derive (44)-(48) for moving dipole sources. From the definition of $f_{n,\varepsilon}$ and the assumption that $m_{k,z}(t)\equiv 0$, we have

[TABLE]

Since $\eta_{\varepsilon}$ is the standard mollifier function, we establish (44).

For ${\mathcal{F}}(g_{n})(\tau),\ {\mathcal{F}}(h_{n})(\tau),\ {\mathcal{F}}(i_{n})(\tau),\ {\mathcal{F}}(j_{n})(\tau)$, changing variable from $t$ to $s$ and integration by parts, we obtain the following form:

[TABLE]

[TABLE]

Then, using similar computation as (54), we establish (45)-(48).

Appendix B. A proof of equation (42)

In this appendix, we set $K(\tau)=K$ for simplicity, and omit the argument $(t_{j}(\tau))$ in $\boldsymbol{p}_{j}(t_{j}(\tau))$.

In $\det\tilde{V}(\tau)$, subtract 1st column from 2nd column, and $(K+1)$-th column from $(K+2)$-th column, we obtain

[TABLE]

Since the 2nd and $(K+2)$-th columns of (55) have the factor $(p_{2,xy}-p_{1,xy})$, we can factorise $\det\tilde{V}(\tau)$ as

[TABLE]

where the $k$-th components of 2$K$-vectors ${\bm\sigma}_{2}$ and ${\bm\sigma}_{2}^{\prime}$ are expressed by

[TABLE]

Also we can derive factor $(p_{2,xy}-p_{1,xy})$ two times by the following two steps

• •

Subtract $(K+1)$-th column from the 2nd column, then 2nd column has the factor $(p_{2,xy}-p_{1,xy})$. Hence, we can factorise by $(p_{2,xy}-p_{1,xy})$.

• •

After above step, subtract twice of 2nd column from $(K+2)$-th column, then $(K+2)$-th column has the factor $(p_{2,xy}-p_{1,xy})$.

Therefore, $\det\tilde{V}(\tau)$ has the factor $(p_{2,xy}-p_{1,xy})^{4}$.

By similar computation, we can derive the factor $(p_{j,xy}-p_{k,xy})^{4}$ for any $j>k$. It is easily see that the largest order of $p_{k,xy}$ in $\det\tilde{V}(\tau)$ is $4(K-1)$ for any $k$, then we can factorise $\det\tilde{V}(\tau)$ as

[TABLE]

where $C_{K}$ is a constant which depends only on $K$.

Let us determine the constant $C_{K}$. By replacing columns $(K-1)$ times and rows $2(2K-2)$ times in $\det\tilde{V}(\tau)$, we obtain

[TABLE]

where

[TABLE]

Therefore the coefficient of the term $(p_{1,xy})^{4K-4}$ is given by

[TABLE]

Applying equation (58) inductively, we obtain

[TABLE]

Here,

[TABLE]

we have $C_{2}=-1$. Then $C_{K}=(-1)^{K(K-1)/2}$, and we complete the proof of (42).