Carleman contraction mapping for a 1D inverse scattering problem with experimental time-dependent data

09/23/2021 ∙ by Thuy T. Le, et al. ∙ UNC Charlotte 0

It is shown that the contraction mapping principle with the involvement of a Carleman Weight Function works for a Coefficient Inverse Problem for a 1D hyperbolic equation. Using a Carleman estimate, the global convergence of the corresponding numerical method is established. Numerical studies for both computationally simulated and experimentally collected data are presented. The experimental part is concerned with the problem of computing dielectric constants of explosive-like targets in the standoff mode using severely underdetermined data.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

We consider in this paper a 1D Coefficient Inverse Problem (CIP) for a hyperbolic PDE and construct a globally convergent numerical method for it. Unlike all previous works of this group on the convexification method for CIPs, which are cited below in this section, we construct in this paper a sequence of linearized boundary value problems with overdetermined boundary data. Each of these problems is solved using a weighted version of the Quasi-Reversibility Method (QRM). The weight is the Carleman Weight Function (CWF), i.e. this is the function, which is involved as the weight in the Carleman estimate for the corresponding PDE operator, see, e.g. [7, 4, 5, 22, 18, 31, 38, 53] for Carleman estimates. Thus, we call this method “Carleman Quasi-Reversibility Method” (CQRM). The key convergence estimate for this sequence is similar with the one of the conventional contraction mapping principle, see the first item of Remarks 8.1 in section 8. This explains the title of our paper. That estimate implies the global convergence of that sequence to the true solution of our CIP. Furthermore, that estimate implies a rapid convergence. As a result, computations are done here in real time, also, see Remark 11.1 in section 11 as well as section 12. The real time computation is obviously important for our target application, which is described below. The numerical method of this paper can be considered as the second generation of the convexification method for CIPs.

The convexification method ends up with the minimization of a weighted globally strictly convex Tikhonov functional, in which the weight is the CWF. The convexification has been developed by our research group since the first two publications in 1995 [20] and 1997 [21]. Only the theory  was presented in those two originating papers. The start for numerical studies of the convexification was given in the paper [1]. Thus, in follow up publications [24, 25, 30, 16, 14, 15, 24, 25, 27, 28, 29, 32, 33, 30, 40, 50, 51] the theoritical results for the convexification are combined with the numerical ones. In particular, some of these papers work with experimentally collected backscattering data [14, 15, 26, 32, 33]. We also refer to the recently published book of Klibanov and Li [31], in which many of these results are described.

The CIP of our paper has a direct application in the problem of the standoff detection and identification of antipersonnel land mines and improvised explosive devices (IEDs). Thus, in the numerical part of this publication, we present results of the numerical performance of our method for both computationally simulated and experimentally collected data for targets mimicking antipersonnel land mines and IEDs. The experimental data were collected in the field, which is more challenging than the case of the data collected in a laboratory [14, 15, 26, 33]. The data used in this paper were collected by the forward looking radar of the US Army Research Laboratory [46]. Since these data were described in the previous publications of our research group [13, 23, 24, 25, 30, 35, 36, 51], then we do not describe them here. Those previous publications were treating these experimental data by several different numerical methods for 1D CIPs. We also note that the CIP of our paper was applied in [32, 33] to the nonlinear synthetic aperture (SAR) imaging, including the cases of experimentally collected data.

Our goal is to compute approximate values of dielectric constants of targets. We point out that our experimental data are severely underdetermined

. Indeed, any target is a 3D object. On the other hand, we have only one experimentally measured time resolved curve for each target. Therefore, we can compute only a sort of an average value of the dielectric constant of each target. This explains the reason of our mathematical modeling here by a 1D hyperbolic PDE rather than by its 3D analog. This mathematical model, along with its frequency domain analog, was considered in all above cited previous publications of our research group, which were working with the experimental data of this paper by a variety of globally convergent inversion techniques.

An estimate of the dielectric constant of a target cannot differentiate with a good certainty between an explosive and a non explosive. However, we hope that these estimates, combined with current targets classification procedures, might potentially result in a decrease of the false alarm rate for those targets. We believe therefore that our results for experimental data might help to address an important application to the standoff detection and identification of antipersonnel land mines and IEDs.

Given a numerical method for a CIP, we call this method globally convergent if:

  1. There is a theorem claiming that this method delivers at least one point in a sufficiently small neighborhood of the true solution of that CIP without relying on a good initial guess about this solution.

  2. This theorem is confirmed numerically.

CIPs are both highly nonlinear and ill-posed. As a result, conventional least squares cost functionals for them are non convex and typically have many local minima and ravines, see, e.g. [49] for a numerical example of this phenomenon. Conventional numerical methods for CIPs are based on minimizations of those functionals, see, e.g. [8, 11, 12]. The goal of the convexification is to avoid local minima and ravines and to provide accurate and reliable solutions of CIPs.

In the convexification, a CWF is a part of a certain least squares cost functional , where is the parameter of the CWF. The presence of the CWF ensures that, with a proper choice of the parameter the functional is strictly convex on a certain bounded convex set of an arbitrary diameter in a Hilbert space. It was proven later in [1] that the existence and uniqueness of the minimizer of on are also guaranteed. Furthermore, the gradient projection method of the minimization of on converges globally to that minimizer if starting at an arbitrary point of . In addition, if the level of noise in the data tends to zero, then the gradient projection method delivers a sequence, which converges to the true solution. Thus, since is an arbitrary number, then this is the global convergence, as long as numerical studies confirm this theory. Furthermore, it was recently established in [33, 40] that the gradient descent method of the minimization of on also converges globally to the true solution.

First, we come up with the same boundary value problem (BVP) for a nonlinear and nonlocal hyperbolic PDE as we did in the convexification method for the same 1D CIP in [30]. The solution of this BVP immediately provides the target unknown coefficient. However, when solving this BVP, rather than constructing a globally strictly convex cost functional for this BVP, as it was done in [30] and other works on the convexification, we construct a sequence of linearized BVPs. Just as the original BVP, each BVP of this sequence involves a hyperbolic operator with non local terms and overdetermined boundary conditions.

Because of non local terms and the overdetermined boundary conditions, we solve each BVP of this sequence by the new version of the QRM, the Carleman Quasi-Reversibility Method (CQRM). Indeed, the classical QRM solves ill-posed/overdetermined BVPs for linear PDEs, see [37] for the originating work as well as, e.g. [19, 31] for updates. Convergence of the QRM-solution to the true solution is proven via an appropriate Carleman estimates [19, 31]. However, a CWF is not used in the optimizing quadratic functional of the conventional QRM of [19, 31]. Unlike this, a CWF is involved in CQRM. It is exactly the presence of the CWF, which enables us to derive the above mentioned convergence estimate, which is an analog of the one of the contraction principle. This allows us, in turn to establish convergence rate of the sequence of CQRM-solutions to the true solution of that BVP. The latter result ensures that our method is a globally convergent one.

Various versions of the second generation of the convexification method involving CQRM were previously published in [2, 3, 45, 44, 39]. In these publications, globally convergent numerical methods for some nonlinear inverse problems are presented. In [2, 3] CIPs for hyperbolic PDEs in are considered. The main difference between [2, 3] and our work is that it is assumed in [2, 3] that one of initial conditions is not vanishing everywhere in the closed domain of interest. In other words, papers [2, 3] work exactly in the framework of the Bukhgeim-Klibanov method, see [6] for the originating work on this method and, e.g. [4, 5, 17, 22, 18, 31, 53] for some follow up publications. On the other hand, in all above cited works on the convexification, including this one, either the initial condition is the function, or a similar requirement holds for the Helmholtz equation. The only exception is the paper [29], which also works within the framework of the method of [6]. We refer to papers [10, 34] for some numerical methods for CIPs with the Dirichlet-to-Neumann map data. In these works the number of free variables in the data exceeds the number of free variables in the unknown coefficient, In our paper Also, in all other above cited works on the convexification.

There is a classical Gelfand-Levitan method [41] for solutions of 1D CIPs for some hyperbolic PDEs. This method does not rely on the optimization, and, therefore, avoids the phenomenon of local minima. It reduces the original CIP to a linear integral equation of the second kind. This equation is called “Gelfand-Levitan equation” (GL). However, the questions about uniqueness and stability of the solution of GL for the case of noisy data is open, see, e.g. Lemma 2.4 in the book [48, Chapter 2]. This lemma is valid only in the case of noiseless data. However, the realistic data are always noisy. In addition, it was demonstrated numerically in [13] that GL cannot work well for exactly the same experimental data as the ones we use in the current paper. On the other hand, it was shown in [24, 25, 51] that the convexification works well with these data. The same is true for the method of the current paper.

Uniqueness and Lipschitz stability theorems of the CIP considered here are well known. Indeed, it was shown in, e.g. [51] that, using the same change of variables as the one considered in section 2 below, one can reduce our CIP to a similar CIP for the equation with the unknown coefficient We refer to [48, Theorem 2.6 of Chapter 2] for the Lipschitz stability estimate for the latter CIP. In addition, both uniqueness and Lipschitz stability results for our CIP actually follow from Theorem 8.1 below as well as from the convergence analysis in two other works of our research group on the convexification method for this CIP [30, 51].

This paper is arranged as follows. In section 2 we state both forward and inverse problems. In section 3 we derive a nonlinear BVP with non local terms. In section 4 we describe iterations of our method to solve that BVP. In section 5 we formulate the Carleman estimate for the principal part of the PDE operator of that BVP. In section 6 we prove the strong convexity of a functional of section 5 on an appropriate bounded set. In section 7 we formulate two methods for finding the unique minimizer of that functional: gradient descent method and gradient projection method and prove the global convergence to the minimizer for each of them. In section 8 we establish the contraction mapping property and prove the global convergence of the method of section 4. In section 9 we formulate two more global convergence theorems, which follow from results of sections 7 and 8. Numerical results with simulated and experimental data are presented in Section 10 and 11 respectively. Concluding remarks are given in section 12.

2 Statements of Forward and Inverse Problems

Below all functions are real valued ones. Let be a known number, be the spatial variable and the function , represents the spatially distributed dielectric constant. We assume that


where is a small number. Let be a positive number. We consider the following Cauchy problem for a 1D hyperbolic PDE with a variable coefficient in the principal part of the operator:


The problem of finding the function from conditions (2.3), (2.4) is our Forward Problem.

Let be the travel time needed for the wave to travel from the point source at to the point ,


By (2.5) the following 1D analog of the eikonal equation is valid:


Let be the Heaviside function,

Lemma 2.1 [30, Lemma 2.1]. For , the function  has the form


where the function  Furthermore,


Lemma 2.2 [Absorbing boundary condition] [30, 50]. Let  be the number in (2.1). Let  and  be two arbitrary numbers. Then the solution  of Forward Problem (2.3), (2.4) satisfies the absorbing boundary condition at  and , i.e.


We are interested in the following inverse problem:

Coefficient Inverse Problem (CIP). Suppose that the following two functions ,  are known:


Determine the function  for , assuming that the number  in (2.1) is known.

Remark 2.1. Note that only the function  can be measured. As to the function  it follows from (2.10) that


We differentiate noisy function using the Tikhonov regularization method [52]. Since this method is well known, we do not describe it here.

3 A Boundary Value Problem for a Nonlinear PDE With Non-Local Terms

We now introduce a change of variable


We will consider the function only for Using (2.6) and (3.1), we obtain


Furthermore, it follows from (2.8)


By (2.5) and (3.3)


Substituting (2.6), (3.3) and (3.4) in (3.2), we obtain


Equation, (3.5) is a nonlinear PDE with respect to the function with nonlocal terms and . We now need to obtain boundary conditions for the function By (2.2) and (2.5) for Hence, (2.11), (2.12) and (3.1) lead to


We will solve equation (3.5) in the rectangle


By (3.1)


By (2.1), (2.2) and (2.5) . Hence, using (2.9) and (3.8), we obtain


It follows from (2.1) and (3.3) that


In addition, we need and we also need to bound the norm from the above. Let be an arbitrary number. We define the set as


We assume below that


By embedding theorem and


where the constant depends only on the domain Using (3.10), we define the function as:


Then the function is piecewise continuously differentiable in and by (3.13) and (3.14)


Thus, (3.5), (3.6), (3.9), (3.14) and (3.16) result in the following BVP for a nonlinear PDE with non-local terms:


Thus, we focus below on the numerical solution of the following problem:

Problem 3.1. Find a function  satisfying conditions (3.17), (3.18), where the function is defined in (3.14).

Suppose that we have solved Problem 3.1. Then, using (3.3) and (3.14), we set


4 Numerical Method for Problem 3.1

4.1 The function

We now find the first approximation for the function Using (2.2), we choose as the first guess for the function Hence, by (3.3),


We now need to find the function To do this, drop the nonlinear third term in the left hand side of equation (3.17) and, using (4.1) and (3.14), set . Then (3.17), (3.18) become:


BVP (4.2), (4.3) has overdetermined boundary conditions (4.3). Typically, QRM works well for BVPs for PDEs with overdetermined boundary conditions [19, 31]. Therefore, we solve BVP (4.2), (4.3) via CQRM. This means that we consider the following minimization problem:

Minimization Problem Number . Assuming (3.12), minimize the functional  on the set


where  is the Carleman Weight Function for the operator [50, 51]


 where  is the parameter, and is the regularization parameter. Both parameters will be chosen later.

Theorem 6.1 guarantees that, for appropriate values of parameters there exists unique minimizer   of the functional

4.2 The function for

Assume that functionals are defined for and their minimizers functions are constructed already, all for the same values of parameters . Replace in (3.17) with with with and with Then problem (3.17), (3.18) becomes a linear one with respect to the function


To solve problem (4.6), (4.7), we consider the following minimization problem:

Minimization Problem Number . Assuming (3.12), minimize the functional  on the set


Suppose that there exists unique minimizer of the functional Then, following (3.19), (3.14) and (3.19), we set


The rest of the analytical part of this paper is devoted to the convergence analysis of the iterative numerical method presented in this section.

5 The Carleman Estimate

In this section we formulate the Carleman estimate, which is the main tool of our construction. This estimate follows from Theorem 3.1 of [30] as well as from (3.15). Let be an arbitrary function and let the function be constructed from the function as in (3.14). Consider the operator

Theorem 5.1 (Carleman estimate [30]). There exists a number  depending only on such that for any  there exists a sufficiently large number  depending only on , such that for all  and for all functions  the following Carleman estimate holds:



Remark 5.1. Here and everywhere below  denotes different constants depending only on listed parameters.

6 Strict Convexity of Functional (4.8) on the Set , Existence and Uniqueness of Its Minimizer

Functional (4.8) is quadratic. We prove in this section that it is strictly convex on the set In addition, we prove existence and uniqueness of its minimizer on this set. Although similar results were proven in many of the above cited publications on the convexification, see, e.g. [30] for the closest one, there are some peculiarities here, which are important for our convergence analysis, see Remarks 6.1, 6.2 below.

Introduce the subspace as:


Denote the scalar product in the space Also, denote


Theorem 6.1. 1. For any set of parameters  and for any the functional  has the Fréchet derivative  The formula for is:


for all

2. This derivative is Lipschitz continuous in , i.e. there exists a constant  such that


3. Let  and be the numbers of Theorem 5.1 Then there exists a sufficiently large constant


depending only on listed parameters such that for all and for all  the functional  is strictly convex on the set  More precisely, let be an arbitrary function and also let the function  Then the following inequality holds:


4. For any there exists unique minimizer


 of the functional  on the set  and the following inequality holds:


where  denotes the scalar product in

Remark 6.1. Even though the expression in the right hand side of (6.4) is linear with respect to the function  we cannot use Riesz theorem here to prove existence and uniqueness of the minimizer  at which  Rather, all what we can prove is (6.9). This is because we need to ensure that the function .

Proof of Theorem 6.1. Since both function satisfy the same boundary conditions, then


By (4.8) and (6.10)


The expression in the first line of (6.11) coincides with expression in the right hand side of (6.4). In fact, this is a bounded linear functional mapping in . Therefore, by Riesz theorem, there exists unique function such that

In addition, it is clear from (6.11) and (LABEL:5.9) that



is the Fréchet derivative of the functional at the point and the right hand side of (6.4) indeed represents Estimate (6.5) obviously follows from (6.4).

We now prove strict convexity estimate (6.7). To do this, we apply Carleman estimate (5.1) to the third line of (6.11). We obtain


By trace theorem there exists a constant depending only on the domain such that

Since the regularization parameter then we can choose so large that

Hence, for these values of the expression in the last line of (6.13) can be estimated from the below as:


Hence, (6.13) and (6.14) imply

This proves (6.7). The existence and uniqueness of the minimizer and inequality (6.9) follow from (6.7) as well as from a combination of Lemma 2.1 and Theorem 2.1 of [1], also see [42, Chapter 10, section 3].

Remark 6.2. Since the functional  is quadratic, then its strict convexity on the whole space  follows immediately from the presence of the regularization term  in it. However, in addition to the claim of its strict convexity, we actually need in our convergence analysis below those terms in the right hand side of the strict convexity estimate (6.7), which are different from the term  These terms are provided by Carleman estimate (6.7). The condition  of Theorem 6.1 is imposed to dominate the negative term in the last line of (6.7).

7 How to Find the Minimizer

Since we search for the minimizer of functional (4.8) on the bounded set rather than on the whole space then we cannot just use Riesz theorem to find this minimizer, also see [1] and [42, Chapter 10, section 3] for the case of finding a minimizer of a strictly convex functional on a bounded set. Two ways of finding the minimizer are described in this section.

7.1 Gradient descent method

Keeping in mind (3.12), choose an arbitrary function We arrange the gradient descent method of the minimization of functional (4.8) as follows:


where is a small number, which is chosen later. It is important to note that since functions then boundary conditions (4.7) are kept the same for all functions Also, using (3.14) and (4.9), we set


Theorem 7.1 claims the global convergence of the gradient descent method (7.1) in the case when

Theorem 7.1. Let the number  be the one defined in (6.6). Let . For this value of  let  be the unique minimizer of the functional  on the set  with the regularization parameter  (Theorem 6.1) Assume that the function  For each choose the starting point of the gradient descent method (7.1) as Then, there exists a number  such that for any  all functions  Furthermore, there exists a number  such that the following convergence estimates are valid:


Proof. Estimate (7.4) follows immediately from [33, Theorem 4.6] combined with “corrections” of functions via (3.14). Estimate (7.5) follows immediately from trace theorem, (3.14) and (7.2)-(7.4).

7.2 Gradient projection method

Suppose now that there is no information on whether or not the function   In this case we construct the gradient projection method. We introduce the function since it is easy to construct the projection operator on a ball with the center at

Suppose that there exists a function such that


This function exists, for example, in the case when the function Indeed, as one of examples of this function (among many), one can take, e.g.

where the function is such that

The existence of such functions is well known from the Analysis course. Denote


Then (4.6), (4.7) become: