Discrete regularization and convergence of the inverse problem for 1+1 dimensional wave equation

TL;DR

This paper introduces a discrete regularization method for solving the inverse boundary value problem of the 1+1 dimensional wave equation, enabling the recovery of wave speed from boundary data with quantifiable error bounds.

Contribution

It proposes a novel discrete regularization strategy for the inverse wave problem and establishes convergence with explicit error estimates.

Findings

Achieved a Hölder-type stability estimate for wave speed reconstruction.

Developed a regularization method applicable to noisy boundary data.

Provided theoretical guarantees for the convergence of the approximation.

Abstract

An inverse boundary value problem for the 1+1 dimensional wave equation $(\partial_t^2 - c(x)^2 \partial_x^2)u(x,t)=0,\quad x\in\mathbb{R}_+$ is considered. We give a discrete regularization strategy to recover wave speed $c(x)$ when we are given the boundary value of the wave, $u(0,t)$, that is produced by a single pulse-like source. The regularization strategy gives an approximative wave speed $\widetilde c$, satisfying a H\"older type estimate $\| \widetilde c-c\|\leq C \epsilon^{\gamma}$, where $\epsilon$ is the noise level.

Equations

$$(\frac{{\partial }^2}{\partial t^2}-c(x{)}^2\frac{{\partial }^2}{\partial x^2})u(t,x)=0,$$

$${\Vert c\Vert }_{C^k(M)}=\sum _{p=0}^k{\sup }_{x\in (0,\infty )}|\frac{{\partial }^{p}c}{\partial x^p}(x)|,$$

$$C_b^k(M)=\{c\in C^k(M);{\Vert c\Vert }_{C^k(M)}<\infty \}.$$

$${\mathcal{V}}^k=$$

$${\Vert c\Vert }_{C^k(M)}\le m,c-1\in C_0^k([L_0,L_1])\}.$$

$$T>\frac{L_1}{C_0}.$$

$$(\frac{{\partial }^2}{\partial t^2}-c(x{)}^2\frac{{\partial }^2}{\partial x^2})u(x,t)=0\text{in M×(0,2T)},$$

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

$$u{|}_{t=0}=0,{\partial }_{t}u{|}_{t=0}=0,$$

$$\Lambda :L^2(0,2T)\to L^2(0,2T),\Lambda f=u^f{|}_{x=0}.$$

$$\mathcal{L}(E)=\{A:E\to E;A\text{ is linear and continuous}\}.$$

$$\mathcal{A}:\mathcal{D}(\mathcal{A})\subset Z\to \mathcal{R}(\mathcal{A})\subset Y,\mathcal{A}(c)={\Lambda }_c.$$

$$\displaystyle\phi_{n,N}(t)=\Big{(}\frac{N}{T}\Big{)}^{\frac{1}{2}}1_{[\frac{(n-1)T}{N},\frac{nT}{N})}(t),\quad t\in[0,2T).$$

$$\displaystyle\mathcal{P}^{N}=span\big{\{}\phi_{1,N},...,\phi_{2N,N}\big{\}}\subset L^{2}(0,2T)$$

$$P^N:L^2(0,2T)\to {\mathcal{P}}^N,P^N(f)=\sum _{j=1}^{2N}⟨f,{\varphi }_{j,N}{⟩}_{L^2(0,2T)}{\varphi }_{j,N}(t).$$

$${\Lambda }_N=P^N\Lambda P^N.$$

$${\mathcal{B}}_E(\mathcal{H},ϵ)=\{\tilde{\mathcal{H}}\in E:{\Vert \tilde{\mathcal{H}}-\mathcal{H}\Vert }_E<ϵ\}.$$

$$\displaystyle l_{0}(\epsilon_{0})=\big{\lfloor}\frac{4}{7}\log_{2}\epsilon_{0}^{-1}\big{\rfloor}$$

$$N_0({ϵ}_0)=2^{l_0}.$$

$$H(t)=\{\begin{array}{ll}1,& t\ge 0,\\ 0,& t<0.\end{array}$$

$${\tilde{m}}_{N,{ϵ}_0}=P^N\Lambda H+n_{N,{ϵ}_0},$$

$${\mathcal{A}}_0:\mathcal{D}({\mathcal{A}}_0)\subset Z\to \mathcal{R}({\mathcal{A}}_0)\subset L^2(0,2T),{\mathcal{A}}_0(c)=\mathcal{A}(c)H=\Lambda H,$$

$${\alpha }_0({ϵ}_0)=a_0{ϵ}_0^{\frac{4}{45}},$$

$$\displaystyle\mathop{\rm sup}\Big{\{}\left\|\mathcal{R}^{(0)}_{N_{0},\alpha_{0}}\widetilde{m}_{N,\epsilon_{0}}-c\right\|_{Z}:\widetilde{m}_{N,\epsilon_{0}}\in\mathcal{P}^{N},N=2^{l}\geq N_{0}(\epsilon_{0}),$$

$$\displaystyle\qquad\left\|\widetilde{m}_{N,\epsilon_{0}}-P^{N}\Lambda H\right\|_{L^{2}(0,2T)}\leq\epsilon_{0}\Big{\}}\leq C\epsilon_{0}^{\gamma_{0}},$$

$${\mathcal{E}}_{N_1}:{\mathcal{P}}^{N_1}\to {\mathcal{P}}^{N_1},$$

$${\tilde{\Lambda }}_{N_1}:L^2(0,2T)\to L^2(0,2T),$$

$${\tilde{\Lambda }}_{N_1}f=\{\begin{array}{ll}{\Lambda }_{N_1}f+{\mathcal{E}}_{N_1}f,& f\in {\mathcal{P}}^{N_1},\\ 0,& f\in ({\mathcal{P}}^{N_1}{)}^{\perp }.\end{array}$$

$${\alpha }_1({ϵ}_1)=a_1{ϵ}_1^{\frac{4}{9}}$$

$$\displaystyle\mathop{\rm sup}\Big{\{}\left\|\mathcal{R}^{(1)}_{N_{1},\alpha_{1}}\widetilde{\Lambda}_{N_{1}}-c\right\|_{Z}:N_{1}\geq\epsilon_{1}^{-4},\quad\widetilde{\Lambda}_{N_{1}}\in\mathcal{L}(\mathcal{P}^{N_{1}}),$$

$$\displaystyle\qquad\left\|\widetilde{\Lambda}_{N_{1}}-\Lambda_{N_{1}}\right\|_{{Y}}\leq\epsilon_{1}\Big{\}}\leq C\epsilon_{1}^{\gamma_{1}},$$

$$1_E(x)=\{\begin{array}{ll}1,& x\in E,\\ 0,& \text{otherwise}.\end{array}$$

$$J:L^2(0,2T)\to L^2(0,2T),Jf(t)=\frac{1}{2}{\int }_0^{2T}{1}_{▲}(t,s)f(s)ds,$$

Full text

Keywords: Inverse problem, regularization theory, wave equation, discretization.

1. Introduction

We consider an inverse boundary value problem for the wave equation

[TABLE]

and introduce a discrete regularization strategy to recover the sound speed $c(x)$ by using the knowledge of perturbed and discetized Neumann-to-Dirichlet map $\widetilde{\Lambda}_{N_{1}}$. Our approach is based on the Boundary Control method [6, 11, 62].

A variant of the Boundary Control method, called the iterative time-reversal control method, was introduced in [14]. The method was later modified in [20] to focus the energy of a wave at a fixed time and in [53] to solve an inverse obstacle problem for a wave equation. In [40] we introduced a modification to the iterative time-reversal control method that is tailored for the 1+1 dimensional wave equation.

The novelty in this paper is that we analyze the effect of the discretization in the regularized solution of the inverse problem. We give a direct discrete regularization method for the non-linear inverse problem for the wave equation. The result contains an explicit (but not necessarily optimal) convergence rate.

By referring to direct methods for non-linear problems we mean the explicit construction of non-linear map to solve the problem without resorting to a local optimization method. In our case the map is given by (43), shown below. The advantage of direct approaches is that they do not suffer from the possibility that the algorithm converges to a local minimum. In particular, they do not require a priori knowledge that the solution is in a small neighbourhood of a given function.

Classical abstract regularization theory is explained in [22]. The iterative regularization of both linear and non-linear inverse problems and convergence rates are discussed in a Hilbert space setting in [15, 24, 26, 49, 51] and in a Banach space setting in [25, 30, 31, 36, 55, 56, 57]. In section 3.5 we compare our regularization strategy to Morozov’s discrepancy principle (MDP). In the context of abstract regularization theory, this principle has been discussed, e.g., in [58].

There are currently only a few regularized direct methods for non-linear inverse problems. An example is a regularisation algorithm for the inverse problem for the conductivity equation in [39]. Also, a direct regularized inversion for blind deconvolution is presented in [28].

2. Regularization Strategy

2.1. Continuity of forward map

We define

[TABLE]

where $M$ denotes the half axis $M=[0,\infty)\subset\mathbb{R}$. We denote the set of bounded $C^{k}(M)$ -functions by

[TABLE]

Let $C_{0},C_{1},{L_{0},L_{1}},m>0$ and define the space of k times differentiable velocity functions

[TABLE]

Here $C_{0}^{k}([L_{0},L_{1}])$ is the subspace of functions in $C^{k}_{b}(M)$ that are supported on $[L_{0},L_{1}]$. Let

[TABLE]

For $c\in{\mathcal{V}^{2}}$ and $f\in L^{2}(0,2T)$, the boundary value problem

[TABLE]

has a unique solution $u=u^{f}\in H^{1}({M\times(0,2T)})$. Using this solution we define the Neumann-to-Dirichlet operator $\Lambda=\Lambda_{c}$,

[TABLE]

For a Banach space $E$ we define

[TABLE]

Let ${Z=C^{2}_{b}(M)}$ and $Y=\mathcal{L}(L^{2}(0,2T))$. The operator $\mathcal{A}$ is defined in the domain $\mathcal{D}(\mathcal{A})=\mathcal{V}^{3}$ by setting

[TABLE]

The notation in (6) means that the range $\mathcal{R}(\mathcal{A})=\mathcal{A}(\mathcal{V}^{3})$ and the domain $\mathcal{D}(\mathcal{A})$ are equipped with the topologies of $Y$ and $Z$, respectively. Note that the maps (5) and (6) are continuous (see [40]).

2.2. Regularization strategies with discretizatized measurements

Let $T$ be as in (3) and $N\in\mathbb{Z}_{+}$. For $n\in\{1,2,3,...,2N\}$ we define the basis functions as

[TABLE]

Note that the functions ${\phi_{n,N}}$ are orthonormal in $L^{2}(0,2T)$. Having (7) we define the space of piecewise constant functions as

[TABLE]

and an orthogonal projection as

[TABLE]

Let $\Lambda$ be as in (5). Using (9) we define

[TABLE]

Let $E$ be a Banach space and $\mathcal{H}\in E$. We denote

[TABLE]

2.2.1. A model for a single discrete and noisy measurement

Let $\epsilon_{0}>0$ and we define

[TABLE]

and

[TABLE]

Let $N=2^{l}\geq N_{0}$, where $l\in\mathbb{Z}_{+}$. Let $P^{N}$ be as in (9) and $\Lambda$ be as in (5). Let us define

[TABLE]

Let us define

[TABLE]

in which $n_{N,\epsilon_{0}}\in\mathcal{P}^{N}$ represents the error and $\left\|n_{N,\epsilon_{0}}\right\|_{L^{2}(0,2T)}\leq\epsilon_{0}$. We consider the quantity $(\epsilon_{0},N_{0},\widetilde{m}_{N,\epsilon_{0}})$ that we call a measurement. Let $\mathcal{A}$ be as in (6) and $H$ be as in (14). We define

[TABLE]

where $\mathcal{D}(\mathcal{A}_{0})=\mathcal{V}^{3}$ . Our main result on the reconstruction of $c(x)$ from the measurement $(\epsilon_{0},N_{0},\widetilde{m}_{N,\epsilon_{0}})$ is given by the following theorem.

Theorem 1**.**

For the operator $\mathcal{A}_{0}:\mathcal{\mathcal{missing}}V^{3}\subset Z\to L^{2}(0,2T)$, there exists an admissible regularization strategy $\mathcal{R}^{(0)}_{N_{0},\alpha_{0}}$ with the choice of parameter

[TABLE]

satisfying the following: For every $c\in{\mathcal{V}^{3}}$ there are $\widetilde{\epsilon}_{0}>$, $a_{0}>0$, and $C>0$ such that

[TABLE]

for all $\epsilon_{0}\in(0,\widetilde{\epsilon}_{0})$. Here $\gamma_{0}=\frac{1}{270}$ and $N_{0}(\epsilon_{0})$ is as in (13).

An explicit bound $\widetilde{\epsilon}_{0}$ and the value for constant $a_{0}$ are given in the proof. The proof of Theorem 1 is given in Section 2.5.

2.2.2. A model for several discrete and noisy measurements

Let

[TABLE]

where $N_{1}\in\mathbb{Z}_{+}$. Having $\Lambda_{N_{1}}$ as in (10) we define a discrete and noisy measurement operator

[TABLE]

Note that $\mathcal{P}^{N_{1}}\cup(\mathcal{P}^{N_{1}})^{\perp}=L^{2}(0,2T)$. With data corresponding to several boundary measurements $(\epsilon_{1},N_{1},\widetilde{\Lambda}_{N_{1}})$ we get the following results with improved error estimates.

Theorem 2**.**

For the operator $\mathcal{A}:\mathcal{\mathcal{missing}}V^{3}\subset Z\to Y$, there exists an admissible regularization strategy $\mathcal{R}^{(1)}_{N_{1},\alpha_{1}}$ with the choice of parameter

[TABLE]

that satisfies the following: For every $c\in{\mathcal{V}^{3}}$ there are $\widetilde{\epsilon}_{1}>0$, $a_{1}>0$, and $C>0$ such that

[TABLE]

for all $\epsilon_{1}\in(0,\widetilde{\epsilon}_{1})$. Here $\gamma_{1}=\frac{1}{54}$.

An explicit bound $\widetilde{\epsilon}_{1}$ and the value for constant $a_{1}$ are given in the proof. The proof of Theorem 2 is given in Section 2.5. We will give explicit choices for $\mathcal{R}^{(0)}_{N_{0},\alpha_{0}}$ in formula (84) and for $\mathcal{R}^{(1)}_{N_{1},\alpha_{1}}$ in formula (101) below. For the convenience of the reader we give a short summary of the regularization strategy. Assume that we are given $\widetilde{\Lambda}_{N}\in\mathcal{L}(\mathcal{P}^{N_{1}})\subset Y$, that is, the discrete Neumann-to-Dirichlet map for the unknown wave speed $c(x)$ with measurements errors. Then the regularization strategy is obtained by the following steps:

1. (1)

Using operator $\widetilde{\Lambda}_{N_{1}}$ in (71) we constructed a source that produces a wave such that $u^{f_{\alpha,r}}(t,x)|_{t=T}$ is close to the indicator function $1_{\mathcal{M}(r)}(x)$ of the domain of influence ${\mathcal{M}(r)}$—see (27). 2. (2)

Using sources $\widetilde{f}^{N}_{\alpha,r}$ we approximately compute the volumes $V(r)$ of the domains of influences—see (28), (70), and Proposition 3. 3. (3)

Using finite differences we compute approximate values of the derivatives of the volumes of the domain influences $\partial_{r}V(r)$—see (80). 4. (4)

We interpolate the obtained values of $\partial_{r}V(r)$. This determines the approximate values of the wave speed $v(r)$ in the travel time coordinates—see (29) and (81). 5. (5)

Finally, we change the coordinates from travel time coordinates to Euclidean coordinates to obtain the approximate values of the wave speed $c(x)$ for $x\in M$—see (82).

2.3. Previous literature

From the point of view of uniqueness questions, the inverse problem for the 1+1 dimensional wave equation is equivalent to the one-dimensional inverse boundary spectral problem. The latter problem was thoroughly studied in the 1950s [23, 41, 50] and we refer to [29, pp. 65–67] for a historical overview. In the 1960s Blagoveščenskiĭ [17, 18] developed an approach to solving the inverse problem for the 1+1 dimensional wave equation without reducing the problem to the inverse boundary spectral problem. This and later dynamical methods have the advantage over spectral methods that they only require data on a finite time interval. Applications of one-dimensional inverse probems have been discussed widely in [16, 29, 33].

The method in the present paper is a variant of the Boundary Control method that was pioneered by M. Belishev [6] and developed by M. Belishev and Y. Kurylev [10, 11] in the late 1980s and early 1990s. Of crucial importance for the method is the result by D. Tataru [62] concerning a Holmgren-type uniqueness theorem for non-analytic coefficients. The Boundary Control method for multidimensional inverse problems has been summarized in [7, 33], and considered for 1+1 dimensional scalar problems in [9, 12, 40] and for multidimensional scalar problems in [32, 34, 42, 46, 47]. For systems it has been considered in [43, 44, 45]. Stability results for the method have been considered in [1] and [35], and computational implementations in [5, 8, 21, 27, 54]. An application of the method to blockage detection in water pipes is in preparation [64].

The inverse problem for the wave equation can also be solved by using complex geometrical optics solutions. These solutions were developed in the context of elliptic inverse boundary value problems [61], and in [52] they were employed to solve an inverse boundary spectral problem. Local stability results can be proven using (real) geometrical optics solutions [13, 59, 60], and in [48] a stability result was proved by using ideas from the Boundary Control method together with complex geometrical optics solutions.

There is an important method based on Carleman estimates [19], often called the Bukhgeim-Klibanov method after its founders, that can be used to show stability results requiring only a single instance of boundary values, when the initial data for the wave equation is non-vanishing. We mention the interesting recent computational work [2] that is based on this method, and also another reconstruction method that uses a single measurement [3, 4]. This method is based on a reduction to a non-linear integro-differential equation, and there are several papers on how to solve this equation (or an approximate version of it)—see [37, 38] for recent results including computational implementations.

2.4. Notations

We will define $\mathcal{R}^{(1)}_{N_{1},\alpha_{1}}$ as a discrete version of the regulation strategy given in [40]. For that we recall some notation from [40].

We denote the indicator function of a set $E$ by

[TABLE]

We define

[TABLE]

where \blacktriangle=\{(t,s)\in(0,2T)^{2};\ \text{t+s\leq 2T$and$s>t>0}\}. We define the time reversal operator as

[TABLE]

and the projections as

[TABLE]

Using (19), (20), and (21) we define that

[TABLE]

and

[TABLE]

We define a regularized inversion with cutoff as

[TABLE]

where $\eta_{Y}:Y\times(0,\infty)\to\mathbb{R}$ is, for example, a continuous function that satisfies $\eta_{Y}(L,\alpha)=1$ when $d(L,Y^{+})\leq\alpha/4$ and $\eta_{Y}(L,\alpha)=0$ when $d(L,Y^{+})\geq\alpha/2$. Here

[TABLE]

We denote by $\boldsymbol{Z}_{\alpha}$ the lift of ${Z}_{\alpha}$ to $C([0,T],Y)$, that is, $\boldsymbol{Z}_{\alpha}(L)(r)$ is ${Z}_{\alpha}(L(r))$. Moreover, we define $b(t)=1_{(0,T)}(t)(T-t)$ and

[TABLE]

We define the travel time coordinates by

[TABLE]

and the domain of influence

[TABLE]

The function $\tau$ is strictly increasing and we denote its inverse by $\chi$. Moreover, $V(r)$ denotes the volume of $\mathcal{M}(r)$ with respect to the measure $dV=c^{-2}dx$, where $c$ is the speed of sound in (4). From [40, Eq. (21)] we see that

[TABLE]

Moreover, according to [40, Eq. (19), (20)], the speed of sound in travel time coordinates $v=c\circ\chi$ satisfies

[TABLE]

Thus $c$ can be computed from $V$. We will next recall how the formula (29) is regularized in [40]. For small $h>0$ we consider the partition

[TABLE]

where $N_{h}\in\mathbb{N}$ satisfies $T-h\leq N_{h}h<T$. We define a discretized and regularized approximation of the derivative operator $\partial_{r}$ by

[TABLE]

[TABLE]

Let us have

[TABLE]

We define an inversion with a cutoff that takes into account the a priori bounds in (2) by

[TABLE]

We denote by $\boldsymbol{z}$ the lift of $z$ to $L^{\infty}(0,T)$, that is, $\boldsymbol{z}(f)(t)=z(f(t))$. We define the extension by one

[TABLE]

and set $W=E\circ\boldsymbol{z}$. We define

[TABLE]

Note that having $f>0$ in (35) we have it that $\tilde{\chi}(f)$ is a strictly increasing function. Having $L_{0},L_{1}$, as in (2), we define

[TABLE]

[TABLE]

Having $f>0$ and using (35) and (36) we define

[TABLE]

Let us define $\eta\in C^{\infty}(\mathbb{R})$ by

[TABLE]

where the constant $C>0$ is selected so that $\int_{\mathbb{R}}\eta(x)=1$. For $\nu>0$ we define

[TABLE]

By using convolution we define a smooth approximation to a given function $f\in L^{\infty}(\mathbb{R})$ by setting

[TABLE]

Using (23), (24), (25), (31), (34), (37), and (42) we define the family of operators for the regularization strategy used in [40] by

[TABLE]

where $\nu=C\epsilon^{\frac{1}{54}}$, $h=C\epsilon^{\frac{1}{18}}$ and $\alpha(\epsilon)=2^{\frac{13}{9}}T^{\frac{4}{9}}\epsilon^{\frac{4}{9}}$. Note that in [40] we considered perturbations of the Neumann-to-Dirichlet operator of the form

[TABLE]

where $\mathcal{E}\in Y$ models the measurement error and $\left\|\mathcal{E}\right\|_{Y}\leq\epsilon$. Below we will introduce a discretized version of regularization strategy (43) that takes in discretized measurements. To this end we start with auxiliary lemmas.

2.5. The proofs of the main results

Lemma 1**.**

Let $N\in\mathbb{Z}_{+}$. Let $\Lambda$ be as in (5) and $\Lambda_{N}$ as in (10). Then we have

[TABLE]

Here $C=C(T)>0$ depends on $T$.

Proof.

By (5) and the trace theorem we have

[TABLE]

By (9) we have

[TABLE]

Let $f\in H^{2}(0,2T)$. By (9) and $H^{2}(0,2T)\hookrightarrow C^{1}([0,2T])$ we have

[TABLE]

Thus

[TABLE]

By (48) and having $C_{sob}=\left\|I\right\|_{H^{2}(0,2T)\to C^{1}([0,2T])}$ we get

[TABLE]

Using interpolation theory with (46) and (49) we have

[TABLE]

For self-adjoint operators $P^{N}$ and $\Lambda$ we have

[TABLE]

By (45) and (50) we have

[TABLE]

Thus

[TABLE]

By (51), (52), and (53) we have

[TABLE]

∎

Proposition 1**.**

Let $\epsilon>0$, $N\in\mathbb{Z}_{+}$ and $N\geq\epsilon^{-4}$ . Let $\Lambda$ be as in (5) and $\Lambda_{N}$ be as in (10). Assume that $\widetilde{\Lambda}_{N}\in\mathcal{B}_{Y}(\Lambda_{N},\epsilon)$, then

[TABLE]

Here $C_{2}=C_{2}(T)>0$ depends on $T$.

Proof.

Using Lemma 1 and having $N\geq\epsilon^{-4}$ we get

[TABLE]

∎

Let $J$ be as in (19) and using (9) we define

[TABLE]

Lemma 2**.**

Let $N\in\mathbb{Z}_{+}$. Let $J$ be as in (19) and $J^{N}$ be as in (56). Then we have

[TABLE]

Here $C=C(T)>0$ depends on $T$.

Proof.

By (19) we have

[TABLE]

Using interpolation with (46) and (49) we have

[TABLE]

By (57) and (58) we have

[TABLE]

We have

[TABLE]

By (57), (59), and (60) we get

[TABLE]

∎

Let $J$ be as in (19), $J^{N}$ be as in (56), and $R$ be as in (20). We define

[TABLE]

By (21) and (62) we define

[TABLE]

Proposition 2**.**

Let $\epsilon\in(0,1)$, $\kappa_{3}>0$, and $N\geq\kappa_{3}\epsilon^{-4}$. Let $\Lambda$ be as in (5), $\boldsymbol{H}$ be as in (23), and $\boldsymbol{H}^{N}$ be as in (63). Assume that $\widetilde{\Lambda}_{N}\in\mathcal{B}_{Y}(\Lambda,\epsilon)$, then we have

[TABLE]

Here $C_{3}=C_{3}(T,\kappa_{3})>0$ depends on $T$ and $\kappa_{3}$.

Proof.

By (22) and (62) we get

[TABLE]

Using (20), we have $\left\|R\right\|_{Y}\leq 1$. Using [40, Theorem 5], we have $\left\|\Lambda\right\|_{Y}\leq M_{1}<\infty$. By (64) and Lemma 2 we have

[TABLE]

By [40, Proposition 1] we have

[TABLE]

By (65) and (66), when $\epsilon\in(0,1)$ and $N\geq\kappa_{3}\epsilon^{-4}$ we have

[TABLE]

∎

Let $P_{r}$ be as in (21), $b$ be as in (25), and $P^{N}$ be as in (9). We define

[TABLE]

Lemma 3**.**

Let $L\in C([0,T],Y)$ and $\boldsymbol{S}$ be as defined in (25). Then

[TABLE]

Proof.

By using (25), (68) with the self-adjointness of operator $P^{N}$ we get

[TABLE]

We have $\left\|P_{r}b\right\|_{L^{2}(0,2T)}\leq(\frac{T^{3}}{3})^{\frac{1}{2}}$ and $\left\|P^{N}b-b\right\|_{L^{2}(0,2T)}\leq(\frac{T^{3}}{12N^{2}})^{\frac{1}{2}}.$ ∎

Let $K$ be as in (22) and $K^{N}$ be as in (62). Let $\boldsymbol{H}$ be as in (23) and $\boldsymbol{H}^{N}$ be as in (63). Let $\boldsymbol{Z}_{\alpha}$ be as in (24). For $r\in[0,T]$ and $\alpha>0$ we denote

[TABLE]

Let $\boldsymbol{S}$ be as in (25). Using (69) we denote

[TABLE]

Let $\boldsymbol{S}^{N}$ be as in (68). Using (69) we denote

[TABLE]

Proposition 3**.**

Let $\epsilon\in(0,1)$ and $\kappa_{4}>0$. Let $\alpha=2\epsilon^{\frac{4}{9}}$ and $N\geq\kappa_{4}\epsilon^{-4}$. Let $s_{\alpha}$ be as defined in (70) and $\widetilde{s}^{N}_{\alpha}$ be as defined in (71). Let $\boldsymbol{H}\Lambda(r)\in Y$ be bounded and positive semidefinite for all $r\in[0,T]$. Assume that $\boldsymbol{H}^{N}\widetilde{\Lambda}_{N}\in\mathcal{B}_{C([0,T],Y)}(\boldsymbol{H}\Lambda,\epsilon)$, then we have

[TABLE]

Here $C_{4}=C_{4}(T,\kappa_{3})>0$ depends on $T$ and $\kappa_{3}$.

Proof.

We have

[TABLE]

Let $\boldsymbol{S}$ be as defined in (25). We have $\left\|\boldsymbol{S}\right\|_{C([0,T],Y)\to C([0,T])}\leq\frac{T^{3}}{3}$, see [40, Proposition 2]. Having $\boldsymbol{H}^{N}\widetilde{\Lambda}_{N}\in\mathcal{B}_{C([0,T],Y)}(\boldsymbol{H}\Lambda,\epsilon)$ we get $\left\|\boldsymbol{Z}_{\alpha}(\boldsymbol{H}^{N}\widetilde{\Lambda}_{N})-\boldsymbol{Z}_{\alpha}(\boldsymbol{H}\Lambda)\right\|_{C([0,T],Y)}\leq\frac{1}{2}\epsilon^{\frac{1}{9}}$, see [40, Proposition 2]. Using (70), for the second part of the sum in the right-hand side we get

[TABLE]

where we denote $\Omega=C([0,T],Y)\to C([0,T])$. Using (70) and (71) with Lemma 3, for the first part of the sum in the right-hand side we get

[TABLE]

We have

[TABLE]

where we denote $\Omega=C([0,T],Y)$. Using [40, Eq. (32)] we have

[TABLE]

Having $\epsilon\in(0,1)$, $N\geq\kappa_{4}\epsilon^{-4}$ and $\alpha=2\epsilon^{\frac{4}{9}}$ with use of (74), (75), and (76) we get

[TABLE]

By using (72) (73) and (77) we get the estimate. ∎

Proof of Theorem 2.

Let us consider the measurement ($\epsilon_{1},N_{1},\widetilde{\Lambda}_{N_{1}}$). There is $N_{1}\geq\epsilon_{1}^{-4}$ and $\widetilde{\Lambda}_{N_{1}}$, for which $\left\|\widetilde{\Lambda}_{N_{1}}-\Lambda_{N_{1}}\right\|_{Y}\leq\epsilon_{1}$, and Proposition 1 gives us

[TABLE]

Let $\epsilon_{1}\in(0,C_{2}^{-1})$ and $\epsilon=C_{2}\epsilon_{1}$. We choose $\kappa_{3}=C_{2}^{4}$ and thus $N_{1}\geq\epsilon_{1}^{-4}=C_{2}^{4}\epsilon^{-4}$. Having $\widetilde{\Lambda}_{N_{1}}\in\mathcal{B}_{Y}(\Lambda,\epsilon)$, Proposition 2 gives us

[TABLE]

Let $\epsilon_{1}\in(0,C_{2}^{-1}C_{3}^{-1})$ and $\epsilon=C_{2}C_{3}\epsilon_{1}$. We choose $\kappa_{4}=C_{2}^{4}C_{3}^{4}$ and thus $N_{1}\geq\epsilon_{1}^{-4}=C_{2}^{4}C_{3}^{4}\epsilon^{-4}$. Let $\alpha_{1}=2C_{2}^{\frac{4}{9}}C_{3}^{\frac{4}{9}}\epsilon_{1}^{\frac{4}{9}}$. Having (78), Proposition 3 gives us

[TABLE]

Let $V$ be as defined as in (28) and $D_{h}$ be as defined as in (31). Let $\epsilon_{1}\in(0,\widehat{\epsilon})$, where $\widehat{\epsilon}$ as defined in [40, Proposition 4]. Let $\kappa_{5}>0$ and $h=\kappa_{5}\epsilon_{1}^{\frac{1}{18}}$. Having (79) and using [40, Proposition 4] we get

[TABLE]

where $C_{5}=C_{5}(\kappa_{5})$. Note that we use the parameter $\kappa_{5}$ to control the size of discretization in (31). Let $v$ be as defined in (29) and $W$ be as defined in (34). Let us denote $\widetilde{w}^{N_{1}}_{\alpha}=W(D_{h}(\widetilde{s}^{N_{1}}_{\alpha}))$. Having (80) and using [40, Proposition 5] we get

[TABLE]

Let $\nu=\kappa_{6}\epsilon_{1}^{\frac{1}{54}}$ and $c\in\mathcal{V}^{3}$. Let $\Phi$ be as in (37) and $\eta_{\nu}$ be as in (41). Let us denote $\widetilde{c}^{N_{1}}_{\alpha}=\eta_{\nu}\ast\Phi(\widetilde{w}^{N_{1}}_{\alpha})$. Having (81) and using [40, Proposition 6] we get

[TABLE]

where $C_{7}=C_{7}(\kappa_{6})$. Note that we use the parameter $\kappa_{6}$ to control the support of $\eta_{\nu}$ in (41). We define

[TABLE]

By using (24), (31), (34), (37), (42), (63), and (68) we define

[TABLE]

and have the estimate

[TABLE]

when $\epsilon_{1}\in(0,\widetilde{\epsilon}_{1})$ and $a_{1}=C_{7}C_{6}^{\frac{1}{3}}C_{5}^{\frac{1}{3}}C_{4}^{\frac{1}{6}}C_{3}^{\frac{1}{54}}C_{2}^{\frac{1}{54}}$. ∎

Proof of Theorem 1.

Let ${\phi_{j,N_{0}}}$ be as in (7), where $j\in\{1,2,3,...,2N_{0}\}$. Let $r\in[0,2T)$ and we define

[TABLE]

Using (15) and (85) we define

[TABLE]

As \big{\{}\phi_{1,N_{0}},...,\phi_{2N_{0},N_{0}}\big{\}} span $\mathcal{P}^{N_{0}}$ this defines a linear map

[TABLE]

Using (86) and (87) we define a perturbed and discretizatized Neumann-to-Dirichlet operator

[TABLE]

where $\mathcal{P}^{N_{0}}\oplus(\mathcal{P}^{N_{0}})^{\perp}=L^{2}(0,2T)$. Using (15), (86), (87), and (88) we define

[TABLE]

Let $\epsilon_{0}>0$. Let $l_{0}(\epsilon_{0})$ be as in (12) and $N_{0}(\epsilon_{0})$ be as in (13). Let $N=2^{l}\geq N_{0}$, where $l\in\mathbb{Z}_{+}$. Let $\Lambda_{N_{0}}$ be as in (10) and $\widetilde{\Lambda}_{N_{0}}$ be as in (88). Let $P^{N}$ be as in (9) and $\phi_{1,N}$ be as in (7). Let $\widetilde{m}_{N,\epsilon_{0}}$ be as in (15). Assume that $\widetilde{m}_{N,\epsilon_{0}}\in\mathcal{B}_{L^{2}(0,2T)}(P^{N}\Lambda H,\epsilon_{0})$ and $\left\|n_{N,\epsilon_{0}}\right\|_{L^{2}(0,2T)}\leq\epsilon_{0}$. Then

[TABLE]

Here $C_{P_{1}}=C_{P_{1}}(T)>0$ depends on $T$. Having (7), (14), and (85) we get

[TABLE]

As $N=2^{l}\geq N_{0}=2^{l_{0}}$ we have $P^{N_{0}}P^{N}=P^{N_{0}}$. Using $\mathcal{T}_{r}\Lambda=\Lambda\mathcal{T}_{r}$ and $P^{N_{0}}\mathcal{T}_{\frac{T}{N_{0}}}=\mathcal{T}_{\frac{T}{N_{0}}}P^{N_{0}}$ with (15), (86), and (91) we get

[TABLE]

Using $N_{0}\leq\epsilon_{0}^{-\frac{4}{5}}$, $\left\|\mathcal{T}_{\frac{T}{N_{0}}}\right\|_{L^{2}(0,2T)}\leq 1$ and $\left\|n_{N,\epsilon_{0}}\right\|_{L^{2}(0,2T)}\leq\epsilon_{0}$ we get

[TABLE]

Let $f\in L^{2}(0,2T)$. By (9) we get

[TABLE]

Using (7), (86), (91), and (92) we have

[TABLE]

Thus (93) and (95) gives us

[TABLE]

We have $N_{0}\leq\epsilon_{0}^{-\frac{4}{5}}\leq 2N_{0}$ and this proves (90).

We define

[TABLE]

Using $N_{0}\leq\epsilon_{0}^{-\frac{4}{5}}\leq 2N_{0}$ and (97) we get

[TABLE]

By (98), $\widetilde{\Lambda}_{N_{0}}\in\mathcal{B}_{Y}(\Lambda_{N_{0}},\epsilon_{1})$, and Theorem 2 we get

[TABLE]

Note that in the proof of Theorem 2 we have assumed $N_{1}\geq\epsilon_{1}^{-4}$. After replacing this by (98), the proof is identical, only the constants $\widetilde{C}_{2},\widetilde{C}_{3},\widetilde{C}_{4}$ change—see (55), (67), (76), and (82). Using (97) we get

[TABLE]

We define

[TABLE]

where $\widehat{\epsilon}$ can be specified by using [40, Proposition 4]. By using (24), (31), (34), (37), (42), (89), (63), and (68), we define

[TABLE]

and have the estimate

[TABLE]

when $\epsilon_{0}\in(0,\widetilde{\epsilon}_{0})$ and $a_{0}=C_{7}C_{6}^{\frac{1}{3}}C_{5}^{\frac{1}{3}}\widetilde{C}_{4}^{\frac{1}{6}}\widetilde{C}_{3}^{\frac{1}{54}}\widetilde{C}_{2}^{\frac{1}{54}}C_{P_{1}}^{\frac{1}{54}}$. ∎

3. Numerical examples

In this section we describe a computational implementation of the regularization strategy in Theorem 2. We will also compare this with a heuristic variant of MDP—see (116). We will begin by describing how the data—that is, the noisy discretized Neumann-to-Dirichlet map $\widetilde{\Lambda}_{N_{1}}$—is simulated.

3.1. The simulation of measurement data

We choose $T=0.6$ in all the simulations. We use $k$-Wave [63] to solve the boundary value problem (4) with $f=\phi_{1,N_{1}}\in\mathcal{P}_{N_{1}}$, where $N_{1}=2^{10}$, and denote the solution by $u^{(sim)}$. Recall that

[TABLE]

where $h=T/N_{1}$. In order to simulate $u^{(sim)}$ for $2T$ time units, a fine discretization needs to be used, and we choose a regular mesh with $2^{13}$ spatial and $N_{2}=2^{15}$ temporal cells. Then we define the simulated Neumann-to-Dirichlet map, acting on the first basis function,

[TABLE]

where $t_{j}=\frac{(j-1)2T}{N_{2}}$, $j\in\{1,2,...,N_{2}\}$, are the temporal grid points. The output of k-Wave is, of course, only an approximation of $u^{(sim)}$ but we do not analyze this simulation error and use the same notation for both $u^{(sim)}$ and its approximation.

Our primary object of interest is the following discretized version of the Neumann-to-Dirichlet map

[TABLE]

where $\Lambda_{j}=\langle\Lambda^{(sim)}\phi_{1,N_{1}},\phi_{j,N_{1}}\rangle_{L^{2}(0,2T)}$ and $f_{k}$, $k=1,2,\dots,2N_{1}$, are the coefficients of $f$ on the basis of $\mathcal{P}^{N_{1}}$. Observe that $\Lambda_{N_{1}}^{(d)}\phi_{1,N_{1}}$ is simply the projection of $\Lambda^{(sim)}\phi_{1,N_{1}}$ on $\mathcal{P}^{N_{1}}$, and that $\Lambda_{N_{1}}^{(d)}f$ is then defined by using the fact that the wave equation (4) is invariant with respect to translations in time.

We will now describe how the noise is simulated. Consider

[TABLE]

where $n_{j}\in\mathcal{N}(0,1)$, that is, $n_{j}$ is a normally distributed random variable with zero mean and unit variance. We compute a realization of $n$ by using the randn function of MATLAB, and use the same notation for $n$ and its realization. Let $\epsilon_{0}^{(d)}>0$ and define, analogously to (104), a noisy, discretized version of the Neumann-to-Dirichlet map

[TABLE]

where $\widetilde{\Lambda}_{j}=\Lambda_{j}+\widehat{n}_{j}$ and

[TABLE]

Following the formulation of Theorem 2, rather than using $\epsilon_{0}^{(d)}>0$, we prefer to parametrize the noise level in terms of

[TABLE]

In what follows, we will consider the quantity $(\epsilon^{(d)}_{1},\widetilde{\Lambda}_{N_{1}}^{(d)})$, a simulated analogue of the noisy measurements in Section 2.2.2.

3.2. Implementation of the regularization strategy

For the a priori bounds in (2) we use values $C_{0}=0.01$, $C_{1}=10$, $L_{0}=0.01$, and $L_{1}=0.6$. The crux of the regularization strategy $\mathcal{R}^{(1)}_{N_{1},\alpha_{1}}$ is the computation of the inverse in (24). When starting from the simulated measurement $(\epsilon^{(d)}_{1},\widetilde{\Lambda}_{N_{1}}^{(d)})$, the analogue of (24) is to solve $X_{j}$ in the equation

[TABLE]

Here $P_{r}$ is the projection in (21), and we choose $r_{j}=jh$, $j=1,2,\dots,N_{1}$. The choice of the regularization parameter $\alpha_{1}$ is discussed in detail below.

We use the restarted generalized minimal residual (GMRES) method to solve the system of linear equations (109) and choose six as the maximum number of outer iterations and 10 as the number of inner iterations (restarts). We use the initial guess $f=0$ and the tolerance of the method is set to 1e-12.

After this we simply follow the regularization strategy (84), that is, we get an approximation of $c$ by setting

[TABLE]

where $\widetilde{s}^{N_{1}}_{\alpha_{1}}(r_{j})=\langle X_{j},b\rangle_{L^{2}(0,2T)}$, $j=1,2,\dots,N_{1}$. The scaling of $\nu$ is chosen as follows

[TABLE]

In the numerical computations the parameter $h$ was fixed to be $h=\frac{T}{N_{1}}$, that is, the discrete derivative $D_{h}$ was computed in the grid that is used in (102) to represent the basis functions $\phi_{1,N_{1}}$. Observe that this deviates from the theoretical choice $h=C\epsilon_{1}^{\frac{1}{18}}$ used in (43). We will describe next how the regularization parameter $\alpha_{1}$ is chosen, and then we will study how the error $\widetilde{c}^{N_{1}}_{\alpha_{1}}-c$ behaves as function of $\epsilon^{(d)}_{1}$.

3.3. Calibration of the regularization strategy

Recall that in Theorem 2 the choice of regularization parameter is of the form $\alpha_{1}=C_{reg1}(\epsilon_{1}^{(d)})^{p}$, where $p={\frac{4}{9}}$. In particular, the choice is explicit apart from the constant $C_{reg1}$. In this section we choose $C_{reg1}$ so that it gives a good reconstruction of a particular velocity function $c_{c}$—see Figure 1. Then the same constant is used in all the subsequent computational examples.

In the regularization strategy we consider 10 values for measurement errors, as defined in (108)

[TABLE]

and nine values for the multiplicative constant $C_{reg1}=10^{-j}$, $j=1,3,...,9$. Then we consider the error in the reconstruction as a function of $j$,

[TABLE]

where for each error level, the reconstruction $\widetilde{c}^{N_{1}}_{\alpha_{j,k}}$ is computed by (110). These computations are summarized in Figure 2. We see that the choice $j=4$, that is,

[TABLE]

gives a good reconstruction on all the error levels. In what follows we will systematically use this choice.

3.4. Reconstruction results based on the analysis

We will now consider the reconstruction (110), with the choice of regularization parameter (113), in two test cases. We begin with with a smooth velocity function $c_{s}$(see Figure 3), where reconstructions of two different noise levels are shown.

To study the order of convergence of our reconstruction method, we consider 10 noise levels,

[TABLE]

and simulate noisy measurements with five different realizations of the random vector $n$ in (105) at each noise level. The corresponding reconstruction errors $\left\|\widetilde{c}^{N_{1}}_{\alpha_{1}}-c_{s}\right\|_{L^{2}(M)}$ are summarized in Figure 4. Computations suggest that the order of convergence is 0.40. This is better than $\frac{1}{54}$ in Theorem (2).

We also tested the method with a non-smooth velocity function $c_{p}$(see Figure 5), where reconstructions of two different noise levels are shown. This case is not covered by the above analysis, but the reconstruction method is also robust in this case.

3.5. Reconstruction results based on MDP

Here we use a heuristic version of MDP as a parameter choice rule for $\alpha=\alpha(\epsilon)$. Typically MDP is applied to a Tikhonov regularization of the form

[TABLE]

where $F$ is a model for the measurements, $y$ is the data, and $x^{*}$ plays the role of a selection criterion. In our case, the model $F$ corresponds to $c\mapsto\Lambda_{N_{1}}^{c}\phi_{1,N_{1}}$ and $y=\widetilde{\Lambda}_{N_{1}}^{c}\phi_{1,N_{1}}$ gives the measurement data, but we do not cast the inverse problem as a minimization problem and our regularization method is not of the Tikhonov type. In particular, our method does not depend on the choice of the auxiliary parameter $x^{*}$, which can be viewed as an initial guess, and that chooses a local minimum of the non-linear optimization problem (115). Due to these differences, the existing results on MDP do not apply to our method, and (116) below is only a heuristic analogue of the classical MDP. We refer to [58] for a study of MDP in an abstract context of the form (115), with non-linear $F$.

The heuristic principle that we use is as follows. We fix tuning parameters $h>1$ and small $\delta>0$ and search for a regularization parameter $\alpha$ in such a way that the following consistency condition holds:

[TABLE]

Here $\epsilon>0$ is the noise level, $y=\widetilde{\Lambda}_{N_{1}}^{c}\phi_{1,N_{1}}$ is again the measurement data, and $\Lambda_{N_{1}}^{\widetilde{c}^{N_{1}}_{\alpha}}\phi_{1,N_{1}}$ is the corresponding data computed with the velocity function $\widetilde{c}^{N_{1}}_{\alpha}$, given by the reconstruction method. Observe that (116) is a relaxed version of (1.7) in [58].

We choose $h=1.1$ and $\delta=0.01$ and use a bisection search to find $\alpha$. Our implementation was unable to find $\alpha$ satisfying the constraint (116) for noise levels $\epsilon^{(d,k)}>0.02$. For smaller noise levels, the regularization parameters found using the principle are summarized in Figure 6. We see that, with the above choice of tuning parameters, the heuristic MDP always gives a larger regularization parameter than (113). The reconstructions are consistently worse than those produced by the choice (113).

Acknowledgements. We thank Samuli Siltanen for inspiring discussions on the regularization on inverse problems.

L. Oksanen was partly supported by EPSRC, project EP/P01593X/1. J. Korpela and M. Lassas were supported by the Academy of Finland, projects 263235, 273979, 284715, and 312119.