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

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.


Introduction
We consider an inverse boundary value problem for the wave equation x) = 0, and introduce a discrete regularization strategy to recover the sound speed c(x) by using the knowledge of perturbed and discetized Neumannto-Dirichlet map Λ N 1 . Our approach is based on the Boundary Control method [6,11,62].
A variant of the Boundary Control method, called the iterative timereversal 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 (41), 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 nonlinear 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]. Let C 0 , C 1 , L 0 , L 1 , m > 0 and define the space of k times differentiable velocity functions

Regularization Strategy
Here C k 0 ([L 0 , L 1 ]) is the subspace of functions in C k b (M ) that are supported on [L 0 , L 1 ]. Let For c ∈ V 2 and f ∈ L 2 (0, 2T ), the boundary value problem has a unique solution u = u f ∈ H 1 (M × (0, 2T )). Using this solution we define the Neumann-to-Dirichlet operator Λ = Λ c , For a Banach space E we define The notation in (6) means that the range R(A) = A(V 3 ) and the domain D(A) are equipped with the topologies of Y and Z, respectively. Note that the maps (5) and (6) are continuous (see [40]).

Regularization strategies with discretizatized measurements.
Let T be as in (3) and N ∈ Z + . For n ∈ {1, 2, 3, ..., 2N } we define the basis functions as Note that the functions φ n,N are orthonormal in L 2 (0, 2T ). Having (7) we define the space of piecewise constant functions as and an orthogonal projection as Let Λ be as in (5). Using (9) we define Let E be a Banach space and H ∈ E. We denote 2.2.1. A model for a single discrete and noisy measurement. Let 0 > 0 and we define and Let N = 2 l ≥ N 0 , where l ∈ Z + . Let P N be as in (9) and Λ be as in (5). Let us define Let us define in which n N, 0 ∈ P N represents the error and n N, 0 L 2 (0,2T ) ≤ 0 . We consider the quantity ( 0 , N 0 , m N, 0 ) that we call a measurement. Let A be as in (6) and H be as in (14). We define where D(A 0 ) = V 3 . Our main result on the reconstruction of c(x) from the measurement ( 0 , N 0 , m N, 0 ) is given by the following theorem. satisfying the following: For every c ∈ V 3 there are 0 >, a 0 > 0, and for all 0 ∈ (0, 0 ). Here γ 0 = 1 270 and N 0 ( 0 ) is as in (13). An explicit bound 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 where N 1 ∈ Z + . Having Λ N 1 as in (10) we define a discrete and noisy measurement operator With data corresponding to several boundary measurements ( 1 , N 1 , Λ N 1 ) we get the following results with improved error estimates. 1 that satisfies the following: For every c ∈ V 3 there are 1 > 0, a 1 > 0, and C > 0 such that for all 1 ∈ (0, 1 ). Here γ 1 = 1 54 . An explicit bound 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 R (0) N 0 ,α 0 in formula (82) and for R (1) N 1 ,α 1 in formula (99) below. For the convenience of the reader we give a short summary of the regularization strategy. Assume that we are given Λ N ∈ L(P N 1 ) ⊂ 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) Using operator Λ N 1 in (69) we constructed a source that produces a wave such that u fα,r (t, x)| t=T is close to the indicator function 1 M(r) (x) of the domain of influence M(r)-see (27).   [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 onedimensional 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 nonvanishing. 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 R (1) N 1 ,α 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 We define where = {(t, s) ∈ (0, 2T ) 2 ; t + s ≤ 2T and s > t > 0}. We define the time reversal operator as and the projections as Using (19), (20), and (21) we define that and We define a regularized inversion with cutoff as We define the travel time coordinates by and the domain of influence The function τ is strictly increasing and we denote its inverse by χ. Moreover, V (r) denotes the volume of 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 Moreover, according to [40,Eq. (19), (20)], the speed of sound in travel 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 We define a discretized and regularized approximation of the derivative operator ∂ r by Let us have We define an inversion with a cutoff that takes into account the a priori bounds in (2) by We denote by z the lift of z to L ∞ (0, T ), that is, z(f )(t) = z(f (t)). We define the extension by one and set W = E • z. We definẽ Note that having f > 0 in (35) we have it thatχ(f ) is a strictly increasing function. Having L 0 , L 1 , as in (2), we define Having f > 0 and using (35) and (36) we define where the constant C > 0 is selected so that R η(x) = 1. For ν > 0 we define By using convolution we define a smooth approximation to a given function f ∈ L ∞ (R) by setting Using (23), (24), (25), (31), (34), (37), and (40) we define the family of operators for the regularization strategy used in [40] by . Note that in [40] we considered perturbations of the Neumann-to-Dirichlet operator of the form where E ∈ Y models the measurement error and E Y ≤ . Below we will introduce a discretized version of regularization strategy (41) 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 ∈ Z + . Let Λ be as in (5) and Λ N as in (10). Then we have Here C = C(T ) > 0 depends on T .
Proof. Using Lemma 1 and having N ≥ −4 we get Let J be as in (19) and using (9) we define Lemma 2. Let N ∈ Z + . Let J be as in (19) and J N be as in (54). Then we have Here C = C(T ) > 0 depends on T .
Proof. We have Let S be as defined in (25).
where we denote Ω = C([0, T ], Y ) → C([0, T ]). Using (68) and (69) with Lemma 3, for the first part of the sum in the right-hand side we get We have

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 (114). We will begin by describing how the data-that is, the noisy discretized Neumann-to-Dirichlet map Λ 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 = φ 1,N 1 ∈ P N 1 , where N 1 = 2 10 , and denote the solution by u (sim) . Recall that 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, where t j = (j−1)2T N 2 , j ∈ {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 Λ (d) where Λ j = Λ (sim) φ 1,N 1 , φ j,N 1 L 2 (0,2T ) and f k , k = 1, 2, . . . , 2N 1 , are the coefficients of f on the basis of P N 1 . Observe that Λ (d) N 1 φ 1,N 1 is simply the projection of Λ (sim) φ 1,N 1 on P N 1 , and that Λ (d) N 1 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 where n j ∈ 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 where Λ j = Λ j + n j and Following the formulation of Theorem 2, rather than using (d) 0 > 0, we prefer to parametrize the noise level in terms of (d) In what follows, we will consider the quantity ( 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 R (1) N 1 ,α 1 is the computation of the inverse in (24). When starting from the simulated measurement ( , the analogue of (24) is to solve X j in the equation Here P r is the projection in (21), and we choose r j = jh, j = 1, 2, . . . , N 1 . The choice of the regularization parameter α 1 is discussed in detail below.
We use the restarted generalized minimal residual (GMRES) method to solve the system of linear equations (107) 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 (82), that is, we get an approximation of c by setting where s N 1 α 1 (r j ) = X j , b L 2 (0,2T ) , j = 1, 2, . . . , N 1 . The scaling of ν is chosen as follows In the numerical computations the parameter h was fixed to be h = T N 1 , that is, the discrete derivative D h was computed in the grid that is used in (100) to represent the basis functions φ 1,N 1 . Observe that this deviates from the theoretical choice h = C 1 18 1 used in (41). We will describe next how the regularization parameter α 1 is chosen, and then we will study how the error c N 1 α 1 − c behaves as function of (d) Figure 1. The velocity function c c used in the calibration of the regularization strategy.
3.3. Calibration of the regularization strategy. Recall that in Theorem 2 the choice of regularization parameter is of the form 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 (106) (d) 1,k ∈ {k · 10 −2 |k = 1, 2, 3, ..., 10}, 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, where for each error level, the reconstruction c N 1 α j,k is computed by (108). These computations are summarized in Figure 2. We see that the choice j = 4, that is, gives a good reconstruction on all the error levels. In what follows we will systematically use this choice. We also observe that the reconstruction error becomes more sensitive to the choice of C reg1 as the noise level grows.

3.4.
Reconstruction results based on the analysis. We will now consider the reconstruction (108), with the choice of regularization parameter (111), 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, and simulate noisy measurements with five different realizations of the random vector n in (103) at each noise level. The corresponding reconstruction errors c N 1 α 1 − c s L 2 (M ) are summarized in Figure 4. Computations suggest that the order of convergence is 0.40. This is better than 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 α = α( ). Typically MDP is applied to a Tikhonov regularization of the form 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 → Λ c N 1 φ 1,N 1 and y = Λ c N 1 φ 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 (113). Due to these differences, the existing results on MDP do not apply to our method, and (114) 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 (113), with non-linear F .
The heuristic principle that we use is as follows. We fix tuning parameters h > 1 and small δ > 0 and search for a regularization parameter α in such a way that the following consistency condition holds: Here > 0 is the noise level, y = Λ c N 1 φ 1,N 1 is again the measurement data, and Λ c N 1 α N 1 φ 1,N 1 is the corresponding data computed with the velocity function c N 1 α , given by the reconstruction method. Observe that (114) is a relaxed version of (1.7) in [58].
We choose h = 1.1 and δ = 0.01 and use a bisection search to find α. Our implementation was unable to find α satisfying the constraint (114) for noise levels (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 (111). The reconstructions are consistently worse than those produced by the choice (111).