Explicit inversion formulas for the two-dimensional wave equation from Neumann traces

TL;DR

This paper develops explicit inversion formulas to recover initial data of the 2D wave equation from Neumann boundary measurements, providing exact solutions for circular and elliptical domains with numerical validation.

Contribution

It introduces new explicit inversion formulas for the 2D wave equation from Neumann data, including exact solutions for circular and elliptical geometries.

Findings

Derived explicit back-projection inversion formulas.

Achieved exact recovery for circular and elliptical domains.

Numerical results demonstrate accuracy and stability.

Abstract

In this article we study the problem of recovering the initial data of the two-dimensional wave equation from Neumann measurements on a convex domain with smooth boundary in the plane. We derive an explicit inversion formula of a so-called back-projection type and deduce exact inversion formulas for circular and elliptical domains. In addition, for circular domains, we show that the initial data can also be recovered from any linear combination of its solution and its normal derivative on the boundary. Numerical results of our implementation of the derived inversion formulas are presented demonstrating their accuracy and stability.

Equations

$$({\partial }_t^2-\Delta )u(x,t)$$

$$u(x,0)$$

$$({\partial }_{t}u)(x,0)$$

$$au(x,t)+b{\partial }_{\nu }u(x,t),(x,t)\in \partial \Omega \times (0,\infty ),$$

$$b{\partial }_{\nu }u(x,t),(x,t)\in \partial \Omega \times (0,\infty ).$$

$$({\partial }_t^2-\Delta )v(x,t)$$

$$v(x,0)$$

$$({\partial }_{t}v)(x,0)$$

$$({\partial }_t^2-\Delta )v(x,t)$$

$$v(x,0)$$

$$({\partial }_{t}v)(x,0)$$

$$\mathbb{R}\to \mathbb{R}:h\mapsto u(x+h\nu (x),t)$$

$$\mathcal{M}f:{\mathbb{R}}^2\times [0,\infty )\to \mathbb{R}:(x,t)\mapsto \frac{1}{2\pi }{\int }_{{\mathbb{S}}^1}f(x+r\omega )d\sigma (\omega ),$$

$$\mathcal{R}f:{\mathbb{S}}^1\times \mathbb{R}\to \mathbb{R}:(\theta ,s)\mapsto {\int }_{\mathbb{R}}f(s\theta +a{\theta }^{\perp })da,$$

$${\mathcal{H}}_s\phi :{\mathbb{S}}^1\times \mathbb{R}\to \mathbb{R}:(\theta ,s)\mapsto \frac{1}{\pi }\underset{\epsilon \searrow 0}{\lim }{\int }_{(s-\epsilon ,s+\epsilon )}\frac{\phi (\theta ,t)}{s-t}dt$$

$${\mathcal{K}}_{\Omega }f(x):=\frac{1}{8{\pi }^2}{\int }_{\Omega }f(y)\frac{\left({\partial }_s^2{\mathcal{H}}_s\mathcal{R}{\chi }_{\Omega }\right)\left(\tilde{n}(x,y),\tilde{s}(x,y)\right)}{\Vert x-y\Vert }dy,$$

$$u(x,t)$$

$$v(x,t)$$

$${\int }_{{\mathbb{R}}^2}\phi (y)\left({\int }_0^{\infty }({\partial }_r\mathcal{M}f)(y,r)\log \left|r^2-{\Vert x-y\Vert }^2\right|dr\right)dy$$

$$=\frac{1}{2\pi }{\int }_{{\mathbb{R}}^2}\frac{f(y)}{\Vert x-y\Vert }\left({\int }_{\mathbb{R}}\left({\partial }_s\mathcal{R}\phi \right)\left(\tilde{n}(x,y),s\right)\log \left|2\Vert x-y\Vert \left(s-\tilde{s}(x,y)\right)\right|ds\right)dy$$

$$-{\int }_{\Omega }\phi (y)f(y)\log \left({\Vert x-y\Vert }^2\right)dy.$$

$${\int }_0^{\infty }$$

$$=-f(y)\Phi \left(-{\Vert x-y\Vert }^2\right)-2{\int }_0^{\infty }\mathcal{M}f(y,r){\Phi }^{\prime }\left(r^2-{\Vert x-y\Vert }^2\right)rdr.$$

$${\int }_0^{\infty }({\partial }_r\mathcal{M}f)(y,r)\Phi \left(\left|r^2-{\Vert x-y\Vert }^2\right|\right)dr$$

$$=-f(y)\Phi \left(-{\Vert x-y\Vert }^2\right)-\frac{1}{\pi }{\int }_0^{\infty }{\int }_{{\mathbb{S}}^1}f(y+r\omega ,r){\Phi }^{\prime }\left(r^2-{\Vert x-y\Vert }^2\right)rd\sigma (\omega )dr$$

$$=-f(y)\Phi \left(-{\Vert x-y\Vert }^2\right)-\frac{1}{\pi }{\int }_{{\mathbb{R}}^2}f(y+z){\Phi }^{\prime }\left({\Vert z\Vert }^2-{\Vert x-y\Vert }^2\right)dz.$$

$${\int }_{{\mathbb{R}}^2}\phi (y)$$

$$=-{\int }_{{\mathbb{R}}^2}\phi (y)f(y)\Phi \left(-{\Vert x-y\Vert }^2\right)dy-\frac{1}{\pi }{\int }_{\Omega }\phi (y)f(y){\int }_{{\mathbb{R}}^2}{\Phi }^{\prime }\left({\Vert z\Vert }^2-{\Vert x-y\Vert }^2\right)dzdy,$$

$$\frac{1}{2\pi }{\int }_{{\mathbb{R}}^2}\frac{f(y)}{\Vert x-y\Vert }\left({\int }_{\mathbb{R}}\left({\partial }_s\mathcal{R}\phi \right)\left(\tilde{n}(x,y),s\right)\Phi \left(2\Vert x-y\Vert \left(s-\tilde{s}(x,y)\right)\right)ds\right)dy,$$

$${\int }_{\Omega }{\int }_0^{\infty }u(x,t)v(x,t)dtdx=-\frac{1}{8{\pi }^2}{\int }_{\Omega }{\int }_{\Omega }f(x)g(y)\frac{\left({\mathcal{H}}_s\mathcal{R}{\chi }_{\Omega }\right)\left(\tilde{n}(x,y),\tilde{s}(x,y)\right)}{\Vert x-y\Vert }dxdy,$$

$$u(x,t)=f(x)+{\int }_0^t\frac{t}{\sqrt{t^2-r^2}}{\partial }_r\mathcal{M}f(x,r)dr,t\in [0,\infty ).$$

$$\left|u(x,t)\right|$$

$$\left|v(x,t)\right|\le \{\begin{array}{ll}t{\Vert {\partial }_r\mathcal{M}g\Vert }_{\infty }& \text{ if }0<t\le T\\ {\Vert {\partial }_r\mathcal{M}g\Vert }_{\infty }\frac{R^2}{2\sqrt{t^2-R^2}}& \text{ if }t>T.\end{array}$$

$${\int }_0^{\infty }\left|u(x,t)v(x,t)\right|dt={\int }_0^T\left|u(x,t)v(x,t)\right|+{\int }_T^{\infty }\left|u(x,t)v(x,t)\right|\le c_1+c_2{\int }_{T^2-R^2}^{\infty }\frac{1}{t^2}dt<\infty ,$$

$${\int }_{\Omega }{\int }_0^{\infty }u(x,t)v(x,t)dtdx$$

$$=-{\int }_{\Omega }{\int }_0^{\infty }f(x)\log (r_1)r_1\mathcal{M}g(x,r_1)d{r}_{1}dx$$

Full text

Keywords: wave equation, inverse problems, computed tomography, inversion formula, back-projection, Neumann trace

AMS subject classifications: 35R30, 65R32, 65M32

1 Introduction

Inverse source problems for the wave equation form the basis of various practical applications. For example, in photoacoustic tomography (PAT), one is interested to recover the initial pressure distribution from measurements of the induced pressure waves outside of the investigated object. In PAT and other applications, the corresponding initial value problem is given by the wave equation

[TABLE]

where $\Delta$ denotes the Laplacian in the spatial component, $f\colon\mathbb{R}^{d}\to\mathbb{R}$ the initial pressure distribution and $n\geq 2$ the spatial dimension. A typical inverse problem arising in PAT is as follows: Given a bounded domain $\Omega\subset\mathbb{R}^{n}$ with smooth boundary $\partial\Omega$ and the solution of the wave equation (1.1) on $\partial\Omega\times(0,\infty)$, determine the function $f$ in (1.1) where $f$ is assumed to be a smooth function vanishing outside of $\Omega$ [23, 32, 28].

The problem of recovering $f$ from Dirichlet data $\left.u\right|_{\partial\Omega\times(0,\infty)}$ (Dirichlet trace) has been extensively studied in the recent years. Existing techniques for the reconstruction of $f$ are based on Fourier domain algorithms [1, 15, 22, 34], iterative algorithms [12, 4, 17, 3, 6, 30], time reversal [5, 16, 31, 26] and explicit inversion formulas of the back-projection type [8, 10, 24, 11, 19, 25, 2, 27, 29]. In this paper, we focus on the latter class, which provides theoretical inside and serves as basis of efficient numerical algorithms. Inversion formulas for Dirichlet measurements are well-known for special domains. In [8, 7] several exact formulas for spherical domains in odd [8] and even [7] dimension have been derived, whereas in [10, 24, 11] inversion formulas for elliptical domains in the plane [10], 3D space [11] and arbitrary spatial dimension [11] have been derived. In [13, 14] the problem of determining initial data from knowledge of Dirichlet measurements on certain quadratic hypersurfaces, including parabolic surfaces, has also been analyzed. In [20, 21], explicit inversion formulas for the wave equation on the surface of certain polygons, polyhedra and corner-like domains are presented.

Typical ultrasound sensors are often direction dependent and the measured data is a weighted combination of the normal derivative and the pressure [33, 32]. Such measurements can be modeled by

[TABLE]

where $\partial_{\nu}$ denotes the normal derivative of $u\colon\mathbb{R}^{n}\times[0,\infty)\to\mathbb{R}$ along $\partial\Omega$ with respect to the spatial variable and $a,b\neq 0$ are constants. In this paper we study the problem of reconstructing $f$ in (1.1) from data (1.2). The case $b=0$ corresponds to Dirichlet measurements and is the standard case studied in photoacoustic tomography. The case $a=0$ corresponds to Neumann measurements

[TABLE]

To this day, explicit formulas for the inversion of the wave equation from Neumann measurements are hardly known. For mixed measurements (1.2), in [35], series inversion formulas have been derived. In [9], an exact inversion formula in 3D of universal back-projection type for recovering the initial data $g\in C_{c}^{\infty}(\mathbb{B}^{3}_{\rho})$ in

[TABLE]

from Neumann measurements on the boundary of the open origin centered ball $\mathbb{B}^{3}_{\rho}\subset\mathbb{R}^{3}$ with radius $\rho>0$ has been derived.

In the present paper, we establish explicit formulas for the inversion of the two-dimensional wave equation (1.1) from Neumann measurements $\left.\partial_{\nu}u\right|_{\partial\Omega\times(0,\infty)}$ (Neumann trace). In Section 3.1, we derive an explicit formula for the initial data up to a smoothing integral operator for convex domains by the knowledge of (1.3). We will show that the smoothing integral operator vanishes for circular and elliptical domains (see 3.2). In Section 3.3, we derive an inversion formula for circular domains that exactly recovers $f$ from mixed measurements given in (1.2). Numerical results with the derived formulas are presented in Section 4.

2 Notation and auxiliary results

2.1 Notation

Before we present our results, we start with the notation which we are using throughout this article. By $u\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ we denote the solution of the two-dimensional wave equation (1.1) with initial data $(f,0)$ and by $v\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ the solution of the wave equation

[TABLE]

with initial data $(0,g)$, where $f,g\in C_{c}^{\infty}(\Omega)$.

By $\nu(x)\in\mathbb{R}^{2}$ we denote the outward unit normal vector of $\partial\Omega$ in $x\in\partial\Omega$. As already mentioned in the introduction, $\partial_{\nu}u(x,t)$ denotes the normal derivative of $u$ at $(x,t)\in\partial\Omega\times(0,\infty)$ with respect to the spatial variable $x$, that is, the derivative of the function

[TABLE]

at zero. Applying the chain rule in (2.2) we see that $\partial_{\nu}u(x,t)=\left\langle\nabla u(x,t),\nu(x)\right\rangle$ where $\nabla u(x,t)\in\mathbb{R}^{2}$ is the gradient of $u$ with respect to $x$.

For an integrable function $f\colon\mathbb{R}^{2}\to\mathbb{R}$, the spherical mean operator of $f$ is defined as

[TABLE]

where $\mathbb{S}^{1}\coloneqq\left\{x\in\mathbb{R}^{2}\mid\left\lVert x\right\rVert=1\right\}$, with $\left\lVert\,\cdot\,\right\rVert$ being the Euclidian distance is the unit circle, and $\sigma$ the standard surface measure on manifolds. The Radon transform of $f$ is defined in a similar way

[TABLE]

where $\theta^{\perp}$ is a unit vector orthogonal to $\theta$. By

[TABLE]

we denote the Hilbert transform of a function $\varphi\colon\mathbb{S}^{1}\times\mathbb{R}\to\mathbb{R}$ in the second variable. If $\varphi(\theta,\cdot)$ is differentiable for some $\theta\in\mathbb{S}^{1}$, then we also use $\partial_{s}\varphi(\theta,\cdot)$ to denote the derivative of $\varphi$ with respect to the second variable.

The general inversion formula that we derive for an elliptic domain will recover any smooth function $f\in C_{c}^{\infty}(\Omega)$ up to smoothing integral operator

[TABLE]

which coincides with the error operator appearing in [10]. Here $\chi_{\Omega}\colon\mathbb{R}^{2}\to\mathbb{R}$ is the indicator function of $\Omega$ in $\mathbb{R}^{2}$, and for $x,y\in\mathbb{R}^{2}$ with $x\neq y$ we set $\tilde{n}(x,y)\coloneqq(y-x)/\left\lVert x-y\right\rVert$, $\tilde{s}(x,y)\coloneqq(\left\lVert y\right\rVert^{2}-\left\lVert x\right\rVert^{2})/(2\left\lVert x-y\right\rVert)$.

2.2 Auxiliary results

First, we recall the well-known explicit solution formulas (see [18])

[TABLE]

where $f,g\in C_{c}^{\infty}(\Omega)$, $u$ is the solution of (1.1) with initial data $(f,0)$ and $v$ the solution of (2.1) with initial data $(0,g)$.

The first Lemma which we use for derivation of a explicit inversion formulas is an integral identity for the spherical means. It is a corrected version of [10, Lemma 2.1], where the term $\int_{\Omega}\varphi(y)f(y)\log(\left\lVert x-y\right\rVert^{2})dy$ is missing.

Lemma 2.1**.**

Let $f\in C_{c}^{\infty}(\Omega)$ and $\varphi\colon\mathbb{R}^{2}\to\mathbb{R}$ a continuously differentiable function on $\overline{\Omega}$ vanishing outside of $\overline{\Omega}$. Then, for every $x\in\Omega$ the following identity holds

[TABLE]

Proof.

1. (i)

First we show identity (2.6) when $\log\left|\,\cdot\,\right|$ is replaced by a differentiable and integrable function $\Phi\colon\mathbb{R}\to\mathbb{R}$.

Applying integration by parts on the left integral in (2.6) yields

[TABLE]

Using the definition of the spherical mean operator and introducing polar coordinates, we further obtain

[TABLE]

After substituting $z$ with $z-y$ in the last integral, multiplying with $\varphi(y)$ and integrating both sides with respect to $y$, we obtain

[TABLE]

where we used Fubini’s theorem to change the order of the integrals. In [10] it has been shown that second integral on the hand right side equals to

[TABLE]

which shows the desired identity for $\Phi$ in place of $\log\left|\,\cdot\,\right|$. 2. (ii)

Finally, we choose a sequence $(\Phi_{n})_{n\in\mathbb{N}}$ of differentiable and integrable functions converging pointwise to $\log\left|\,\cdot\,\right|$. Thus, applying (i) and Lebesgue’s dominated convergence theorem yield the claimed identity.∎

The next Lemma can be found in [10], where the proof however contains a small error. Below we give a corrected proof based on the corrected Lemma 2.1, which for the convenience of the reader is given in full detail.

Lemma 2.2**.**

Suppose that $f,g\in C_{c}^{\infty}(\Omega)$. Then, we have

[TABLE]

where $u$ is the solution of (1.1) with initial data $(f,0)$ and $v$ the solution of (2.1) with initial data $(0,g)$. Moreover, $\tilde{n}(x,y)\coloneqq(y-x)/\left\lVert x-y\right\rVert$, $\tilde{s}(x,y)\coloneqq(\left\lVert y\right\rVert^{2}-\left\lVert x\right\rVert^{2})/(2\left\lVert x-y\right\rVert)$ for $x\neq y\in\Omega$.

Proof.

1. (i)

First, we show that $u(x,\cdot)v(x,\cdot)$ is integrable on $(0,\infty)$ for any fixed $x\in\Omega$. By applying the Leibnitz-rule and (2.4) we obtain

[TABLE]

Since $\Omega$ is bounded, we can find $R>0$ such that $\partial_{r}\mathcal{M}f(x,t)=\mathcal{M}f(x,t)=0$ for $r\geq R$. Now, let $T>R$ be a fixed positive number. From the above equality we deduce

[TABLE]

for $0<t\leq T$. For $t>T$, we obtain the inequality

[TABLE]

by differentiating under the integral sign in (2.4). Analogously, we deduce from (2.5) the estimates

[TABLE]

Thus, we obtain

[TABLE]

where $c_{1},c_{2}$ are positive constants. This shows the integrability of $u(x,\cdot)v(x,\cdot)$ and in particular the integrability of $uv$ on $\Omega\times(0,\infty)$. 2. (ii)

In the next step we show that

[TABLE]

For that purpose, let $T>0$ be fixed again. From (2.7) we see that

[TABLE]

The right triple-integral can be evaluated to

[TABLE]

where we applied Fubini’s theorem and used the identity

[TABLE]

Next, we use integration by parts and Fubini’s theorem again on (2.9) to obtain

[TABLE]

where the first integral can be extended to

[TABLE]

by a further application of Fubini’s theorem. Thus, we finally have

[TABLE]

Letting $T\to\infty$ and integrating both sides afterwards we see that (2.8) holds.

In the last two steps we reshape both integrals in (2.8) on the right side to prove the final statement. 3. (iii)

Using the definition of spherical mean operator and polar coordinates we observe that second integral on in the right hand side of (2.8) can be evaluated to

[TABLE]

One further application of Fubini’s theorem and substitution of $y$ with $y-x$ lead then to

[TABLE]

and hence, applying Lemma 2.1, we see that this integral coincides with

[TABLE] 4. (iv)

Finally, using polar coordinates and applying the substitution rule for the first integral in (2.8) yields

[TABLE]

Changing the order of integration in the last displayed equation and combining this with Items (ii) and (iii) show the claimed statement.∎

3 Main results

In this section we present and prove our inversion formulas for the inversion of the two-dimensional wave equation from Neumann measurements.

3.1 Formula for Neumann traces on convex domains

The first main result is an explicit reconstruction integral that applies to arbitrary convex domains and yields an exact inversion formula from the Neumann trace up to an explicitly given smoothing integral operator.

Theorem 3.1**.**

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded convex domain with smooth boundary, $f\in C_{c}^{\infty}(\Omega)$ and $u\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ be the solution of the wave equation (1.1) with initial data $(f,0)$ Then, for every $x\in\Omega$, we have

[TABLE]

where the integral operator $\mathcal{K}_{\Omega}$ is defined by (2.3).

We will derive Theorem 3.1 from the following Proposition which is the key ingredient for the derivation of all main results in this paper.

Proposition 3.2**.**

Let $f,g\in C_{c}^{\infty}(\Omega)$. Then the following identity holds:

[TABLE]

where $u$ is the solution of (1.1) with initial data $(f,0)$ and $v$ the solution of (2.1) with initial data $(0,g)$.

Proof.

1. (i)

We first show that

[TABLE]

Application of integration by parts on the two inner integrals yields

[TABLE]

Therefore, the subtraction of both integrals and changing order of integration lead to the desired result. 2. (ii)

In the next step we use Green’s second identity on right inner integral in (3.3) to obtain

[TABLE]

where the right integral can be written as

[TABLE]

Then the product rule yields the relation

[TABLE]

From the divergence theorem we conclude that the second integral on the right side can be evaluated to

[TABLE] 3. (iii)

Now, it remains to show that

[TABLE]

First, we notice that one can easily verify the relation $\Delta uv=v\Delta u+2\left\langle\nabla u,\nabla v\right\rangle+u\Delta v$. Furthermore, from (1.1) and (2.1) we observe that $\Delta u$ and $\Delta v$ are the solutions of (1.1) and (2.1) with respect to $\Delta f$ and $\Delta g$. The same also holds for the gradients $\nabla u$ and $\nabla v$. This observation and Lemma 2.2 imply then

[TABLE]

where

[TABLE]

Next, we apply Green’s first identity and integration by parts on the right integral in (3.4) to deduce

[TABLE]

Thus, we are left to show

[TABLE]

This is straightforward computation: Applying the chain rule on the numerator yields

[TABLE]

Then, from product rule and the relation $(\nabla_{x}+\nabla_{y})\left\lVert x-y\right\rVert^{-1}=0$ we conclude

[TABLE]

One further application of the chain rule and product rule finally yields then

[TABLE]

Now we can proof our main result.

Proof of Theorem 3.1.

First, let $g\in C_{c}^{\infty}(\Omega)$ and $v\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ denote the solution of (2.1) with initial data $(0,g)$. As in Lemma 3.3 stated, we have that

[TABLE]

By using polar coordinates, Equation (2.5) and substitution rule we have

[TABLE]

Thus, inserting this relation into the second integral in (3.5) and applying Fubini’s theorem yields then

[TABLE]

This leads finally to

[TABLE]

Since this identity holds for every test function $g\in C_{c}^{\infty}(\Omega)$, the claimed inversion formula (3.1) holds. ∎

3.2 Exact formula for Neumann traces on ellipses

Next, we present a formula that exactly recovers the initial data $f\colon\mathbb{R}^{2}\to\mathbb{R}$ from Neumann measurements given in (1.3) in the case that $\Omega$ is bounded by an ellipse

[TABLE]

for numbers $e_{1},e_{2}>0$ and an orthogonal transform $Q\in\mathbb{R}^{2\times 2}$.

Theorem 3.3**.**

Let $\Omega\subset\mathbb{R}^{2}$ be a circular or elliptical domain and $f\in C_{c}^{\infty}(\Omega)$. Then, for every $x\in\Omega$ we have

[TABLE]

where $u\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ denotes the solution of the wave equation (1.1) with initial data $(f,0)$.

Proof.

By Theorem 3.1, we are left to show that $\mathcal{K}_{\Omega}(x)=0$ for every $x\in\Omega$. For the proof of this identity, we refer to [10, 11]. ∎

3.3 Exact formula for mixed traces on circular domains

Finally, we show that the initial data $f\colon\mathbb{R}^{2}\to\mathbb{R}$ in (1.1) can be recovered by any linear combination of the solution of the wave equation $u\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ and the normal derivative $\partial_{\nu}u$ on circular domains.

Theorem 3.4**.**

Let $\mathbb{B}_{\rho}^{2}(z)\subset\mathbb{R}^{2}$ be an open ball with radius $\rho>0$ and center $z\in\mathbb{R}^{2}$ in the plane, let $u$ denote the solution of the wave equation (1.1) with initial data $(f,0)$, where $f\in C_{c}^{\infty}(\mathbb{B}_{\rho}^{2}(z))$, and $a\geq 0$ and $b>0$. Then, for every $x\in\mathbb{B}_{\rho}^{2}(z)$, we have,

[TABLE]

The proof of Theorem (3.4) uses the following lemma, which is a range condition for the Dirichlet trace on a circle of the wave equation.

Lemma 3.5**.**

Let $f\in C_{c}^{\infty}(\mathbb{B}_{\rho}^{2}(z))$, where $\rho>0$ and $z\in\mathbb{R}^{2}$, and $u\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ the solution of wave equation (1.1) with initial data $(f,0)$. Then for every $x\in\mathbb{B}_{\rho}^{2}(z)$ we have

[TABLE]

Proof.

It is sufficient to consider case where $\mathbb{B}_{\rho}^{2}(z)=\mathbb{B}^{2}$ is the unit ball centered at the origin. The general case follows from translation and rescaling. Suppose $f\in C_{c}^{\infty}(\mathbb{B}^{2})$ has compact support in the open unit ball $\mathbb{B}^{2}$. For this case, in [7] the identities

[TABLE]

for any $x\in\mathbb{B}^{2}$, have been shown, where $v\colon\mathbb{R}^{2}\times[0,\infty)\to\mathbb{R}$ is the solution of the wave equation (2.1) with initial data $(0,f)$. Thus, by using product rule in (3.9) and identity (3.8) we have for every $x\in\mathbb{B}^{2}$

[TABLE]

which implies

[TABLE]

From (2.4), (2.5) we also have the identity $\partial_{t}v=u$. This shows the claimed identity. ∎

We are now ready to proof the inversion formula for mixed trace as a corollary of the inversion formula for the Neumann trace in Theorem 3.3 and the range condition for the Dirichlet trace derived in Lemma 3.5.

Proof of Theorem 3.4.

The following holds:

• •

From Lemma 3.5 we have for every $x\in\mathbb{B}_{\rho}^{2}(z)$

[TABLE]

• •

From Theorem 3.3 we obtain for every $x\in\mathbb{B}_{\rho}^{2}(z)$

[TABLE]

Adding the last two displayed equations leads the desired result. ∎

4 Numerical experiments

In this section, we give details of the numerical implementation of the derived exact inversion formula (3.6) and show numerical results. The inversion formula (3.7) with mixed measurements can be implemented in the same way. Throughout this section, we assume a circular domain $\mathbb{B}_{\rho}^{2}(z)$ with radius $\rho>0$ and center $z\in\mathbb{R}^{2}$.

4.1 Discretization of initial data and normal derivative

Suppose that we have given discrete data

[TABLE]

of the initial data $f\in C_{c}^{\infty}(\mathbb{B}_{\rho}^{2}(z))$, where $N$ is the image size, $\Delta x\coloneqq 2\rho/(N-1)$ the step size and $x_{i,j}\coloneqq(-\rho+z_{1}+(i-1)\Delta x,-\rho+z_{2}+(j-1)\Delta x)$. Next, we assume that we have $M\coloneqq\lceil 2\rho\pi/\Delta x\rceil$ detector points located at

[TABLE]

where $\varphi_{k}\coloneqq(k-1)\frac{2\pi}{M-1}$ for $k=1,\ldots,M$. The normal derivative is then discretized by

[TABLE]

where $\mathbf{u}[k,l]$ denotes the solution of wave equation with initial data $(f,0)$ at the point $y_{k}$ and time $t_{l}\coloneqq(l-1)\Delta x$, where $L\coloneqq\lfloor T/dx\rfloor+1$ and $T\coloneqq 16\rho$. We remark that we computed the discrete gradient $\nabla\mathbf{u}[k,l]$ by solving the wave equation on the whole grid (with FFT) and using symmetric finite differences to approximate the partial derivatives.

4.2 Implementation of inversion formula

To implement the inversion formula (3.6)

[TABLE]

for $x\in\mathbb{B}_{\rho}^{2}(z)$ we proceed as follows: First, we approximate for all $(k,l)\in\left\{1,\ldots,M\right\}\times\left\{1,\ldots,L\right\}$ the integral

[TABLE]

by product integration. The last integral can be evaluated to

[TABLE]

and hence we have

[TABLE]

The above inversion formula can then be approximated by

[TABLE]

where we used $\mathrm{interp}(\mathbf{A}[k,\ ],\left\lVert x_{i,j}-y_{k}\right\rVert)$ to denote the interpolated value in $\left\lVert x_{i,j}-y_{k}\right\rVert$ of the array $\mathbf{A}[k,\ ]$. The numerical approximation of formula (3.7) with mixed measurements is denoted by $\mathbf{f_{\mathrm{rec,m}}}$.

Now we present numerical results of the discrete data sets $\mathbf{f_{\mathrm{rec,n}}}\in\mathbb{R}^{N\times N}$ and $\mathbf{f_{\mathrm{rec,m}}}\in\mathbb{R}^{N\times N}$.

4.3 Numerical results

The numerical approximation of the inversion formula (3.6) mentioned above has been implemented in Matlab and was tested on the head phantom $\mathbf{f}$ presented in Figure 1 with image size $N=301$, where the support of $f$ is contained in the open unit ball with $\rho=1$ and $z=0$.

The corresponding Dirichlet and Neumann measurements of the head phantom are shown in the first column in Figure 2. We additionally tested inversion formula (3.7) with the weights $a=1$ and $b=2\Delta x$, where the mixed measurements are also presented in the first column in Figure 2. In the second column the simulated data sets with Gaussian noise (with standard deviation equal to $10\%$ of the maximal value) added are shown.

For comparison reasons, we applied the inversion formula

[TABLE]

which exactly recovers the initial data from Dirichlet measurements (see [10, 11]), on Neumann measurements and vice versa. We denote the corresponding discrete version of (4.1) by $\mathbf{f_{\mathrm{rec,d}}}$ (for implementation details, see $\mathbf{f_{\mathrm{rec,n}}}$). The first column in Figure 3 shows the numerical reconstructions obtained by formula (4.1) and the three different simulated data sets, whereas the second column illustrates the numerical results of the inversion formulas (3.6) and (3.7), respectively. The numerical reconstructions applied to the data sets with $10\%$ Gaussian noise are shown in Figure 4. We also tested our reconstruction method with $20\%$ Gaussian noise added to the data to observe the behaviour of the single reconstructions with a higher noise rate (see Figure 5).

In Figure 3 we observe that both implementations of our derived inversion formulas (3.6) and (3.7) show very good results and approximate the head phantom almost perfectly. As we proved in Lemma 3.5, the numerical approximation of the integral (Figure 3 top, right) is close to zero. We point out that the numerical reconstruction $\mathbf{f_{\mathrm{rec,d}}}$ using mixed measurements is also quite good (Figure 3, middle, left).

Finally, we remark that the numerical reconstructions of the head phantom in Figure 4 show that the structure of the head phantom is still visible although Gaussian noise was added to the data. In Figure 5, we also see that numerical implementation of formula (4.1) with mixed measurements even shows a better result as the numerical implementation of the exact formula (3.7). This leads to assumption that the inversion formula (4.1) is also exact for mixed measurements.

5 Conclusion and Outlook

In this paper, we studied the problem of recovering the initial data of the two dimensional wave equation from Neumann traces. We established an explicit inversion formula for convex domains with smooth boundary up to an explicitly computed smoothing integral operator. This integral operator has been seen to vanish for circular and elliptical domains. We also derived an exact reconstruction formula for recovering the initial data from any linear combination of the solution of wave equation and its normal derivative on circular domains. The numerical results in the last section of this article showed that our implementation of the exact inversion formula leads to very good results even with noisy data.

In future work we intend to investigate the problem of determining the initial data of higher dimensional wave equations from Neumann measurements as well as from mixed measurements. We also want to establish explicit inversion formulas for Neumann data and mixed trace for other special domains such as certain quadratic hypersurfaces or even non-convex domains. Moreover, the derivation of inversion formulas from Neumann measurements on a finite time interval $\partial\Omega\times(0,T)$ instead of $\partial\Omega\times(0,\infty)$ could be also of great interest.