Consider the scattering problem for a penetrable inhomogeneous medium in . Below We assume that the scattering object, which occupies a bounded domain in , is characterized by the spatially distributed dielectric constant , where the function has a compact support. In this paper we are particularly interested in the case of that typically appears in applications of non-destructive testing and explosive detection, see for instance [Bukhg1981, Klib95, Klib97] for a similar assumption. Suppose that the object is illuminated by the downward propagating incident plane wave , where , the propagation direction is fixed, and is the wavenumber. Then there arises the scattered wave, and the total wave which is the sum of the incident wave and the scattered wave is governed by the Helmholtz equation as
The scattered wave satisfies the Sommerfeld radiation condition (2), which guarantees that it behaves like a spherically outgoing wave far away from the scattering object. It is well known that the scattering problem (1)–(2) has a unique solution , see [Colto2013]. Now let and consider
Assume that the scatterer as well as the support of the coefficient are contained in , and that these objects do not intersect with . Let and be positive constants such that . We consider the following inverse problem.
Inverse Problem. Assume that we are given the multi-frequency backscatter Cauchy data
where the total wave is generated by incident plane waves with a fixed propagation direction. Determine the function in (1) for , see also Figure 1 for a schematic diagram of the measurement arrangement in the inverse problem.
Uniqueness theorem for this inverse problem can be currently proven only in the case when the right hand side of equation (1) does not equal to zero in . This can be done by the so-called Bukhgeim-Klibanov method, which was originated in [Bukhg1981]
and is based on applications of Carleman estimates to coefficient inverse problems, see, e.g.[Beili2012, KT, Ksurvey] for this method. In addition, uniqueness of the approximate problem can be proven when the truncated Fourier series for (1) is used for that approximation, see, e.g. Theorem 3.2 in [Khoa2019].
This inverse problem belongs to a wider class of coefficient inverse scattering problems which in general aim to recover information about the coefficient (e.g. its support and/or its values) from the knowledge of the scattered wave generated by a number of incident waves. Inverse scattering problems occur in many applications, including non-destructive testing, explosive detection, medical imaging, radar imaging and geophysical exploration. There is a vast literature about theoretical results and numerical solution to inverse scattering problems, see for instance [Colto2013] and references therein. Due to the interest of this paper, we discuss only some numerical methods. The conventional approach is based on the optimization based methods, see, e.g. [Bakus2004, Chavent, Gonch1, Gonch2, Engl1996]. However, it is well known that these methods may suffer from multiple local minima and ravines and their convergence analysis is also unknown in many situations. An important attempt in overcoming the drawbacks of the optimization based methods is the qualitative approach which aims to compute the geometry of the scattering object or the support of the coefficient . We refer to [Cakon2006, Kirsc2008, Colto2013] and references therein for the development of qualitative methods in solving inverse scattering problems. Although one may be able to avoid local minima or the use of advanced a priori information of the solution, still only geometrical information of the scatterer can be reconstructed with qualitative methods. Furthermore, these methods typically require muti-static data which are sometimes not available in practical applications.
The numerical method proposed in this paper is an extended study from a recent new approach called globally convergent numerical methods (GCNM) for solving coefficient inverse problems. We say that a numerical method for a nonlinear ill-posed problem converges globally if there is a rigorous guarantee that it delivers points in a sufficiently small neighborhood of the exact solution of this problem without any advanced knowledge of this neighborhood. The GCNM typically aims to reconstruct a coefficient in an inverse scattering problem using scattering data either for a single direction of the incident plane wave, or, most recently, for many locations of the point source but at a fixed single frequency [Khoa2019]. An interesting feature of GCNM is that in all cases the data are non over-determined. The latter means that the number of free variables in the data equals the number of free variables in the unknown coefficient, . The main advantage of the GCNM is that any version of it avoids the local minimum problem suffered by optimization based methods. Still, any version of GCNM holds the above indicated global convergence property. We refer to [Beili2012, Kliba2018, Nguye2017, Koles2017, Nguye2018] and references therein for theoretical results as well as numerical and experimental data study of the first type of GCNM.
The method of this paper is inspired by the second type of the GCNM, which is called convexification. The development of the convexification has started in 1995 and 1997 by Klibanov [Klib95, Klib97] and continued since then in [BK, KT, KK]. However, those were mostly analytical works since some obstacles existed at that time on the path to the numerical implementation, although see some numerical results for the one-dimensional case in [KT]. Fortunately, in 2017 the work [Bak] has eliminated those obstacles. This generated a number of more recent publications on the convexification [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019], which contain both a rigorous convergence analysis and numerical results. In particular, publications [Kliba2018a, Kliba2019a] are about the verification of the convexification on experimental data.
The central idea of the convexification is to construct of a globally convex weighted Tikhonov-like functional with the Carleman Weight Function in it. The idea of the use of the Carleman Weight Function is an unexpected consequence of the original idea of the Bukhgeim-Klibanov method [Bukhg1981], which was originally aimed only for proofs of uniqueness theorems for coefficient inverse problems. The final step of the convergence analysis of the convexification consists in the proof of the global convergence of the gradient projection method to the exact solution, as long as the level of noise in the data tends to zero. We also refer to another version of the convexification, which has started in the work [Baud1] and has been continued in [Baud2, BBS, LN]. Carleman Weight Functions are also a crucial element of these works. The main difference between these publications and our method is that it is assumed in [Baud1, Baud2, BBS, LN] that the initial condition in a hyperbolic/parabolic PDE is not vanishing in , which unlike our case of the zero right hand side of equation (1).
As to this present paper, our first step is to eliminate the coefficient from the Helmholtz equation using a change of variables. Next, using a truncated Fourier expansion for a function generated by the total wave field, we approximately reformulate the inverse problem as the Cauchy problem for a system of quasilinear elliptic PDEs. The Cauchy boundary data are as follows: on a part of the boundary both Dirichlet and Neumann boundary data are given and no data are given on the rest of the boundary. We then propose a weighted quasi-reversibility method to solve the problem. Inspired by the concept of the convexification, the cost functional in that weighted quasi-reversibility method contains a Carleman Weight Function. This function plays the decisive role in the numerical performance of the method. A method of gradient descent type is explored to find the global minimizer of the cost functional without using any advanced a priori information about it.
Comparing with the above cited recent works on the convexification, the new features of this work are that firstly our algorithm exploits the new Fourier basis in [Kliba2017] to solve a multi-dimensional inverse problem for the Helmholtz equation with multifrequency data and a single direction of the incident plane wave. The latter is mostly related to [Kliba2018a] in which, however, only the one-dimensional version of the inverse problem has been studied. Using the new Fourier basis from [Kliba2017], the 3D inverse problem for the Helmholtz equation with data generated by a moving source (at a fixed frequency) has been also studied in [Khoa2019]. Secondly, a modification during the iteration of the gradient descent method is applied to help the cost functional converge faster. More precisely, we solve the direct problem to update some functions during the iterations of the gradient descent method. Thirdly, unlike the previous works [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019], the reconstruction algorithm proposed in this paper works without using any data completion process, and the numerical study covers challenging cases of scattering objects of different shapes which are characterized by different values of the dielectric constant. We also want to mention that the implementation of the method uses a full term instead of or terms as in the previous works cited above and does not need any cut-off and averaging procedures during the iteration in the algorithm.
The convergence analysis of the method of this paper will be addressed in an incoming publication. To be more precise, we now roughly (i.e. without some details) specify what kind of theorems will be proven in that publication. Analogs of these theorems for the one-dimensional case can be found in [Kliba2018a]. Those theorems claim:
Existence and uniqueness of the minimizer of the functional on for .
Convergence of the gradient projection method of the minimization of the functional on to the exact solution of that approximate coefficient inverse problem if starting from an arbitrary point of . That convergence takes place as long as the level of noise in the data tends to zero.
Since the radius of the ball is an arbitrary number, then this is the desired global convergence property, as defined above. Note that even though the theory requires the parameter to be sufficiently large, the numerical experience of this and all previous publications about the convexification [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019] shows that actually reasonable values of provide accurate solutions of considered inverse problems.
The paper is structured as follows. The second section is dedicated to the formulation of the inverse problem as an approximate quasilinear elliptic PDE system. The numerical reconstruction method for solving the inverse problem is proposed in Section 3. The implementation and numerical examples of the reconstruction method are presented in Section 4. Finally, Section 5 contains a summary discussion of this work.
2 An approximate elliptic PDE formulation
In this section we reduce our inverse problem to the Cauchy problem for a system of quasilinear elliptic PDEs that we will be studying using a quasi-reversibility approach in the next section. The main ideas for deriving the formulation are using truncated Fourier expansion in and eliminating the coefficient from the scattering problem. Setting we first need the following important Fourier basis of that was introduced in [Kliba2017]
Applying the Gram–Schmidt process towe obtain an orthonormal basis which has the following properties, also, see [Kliba2017]:
i) for all
ii) The matrix , where and
is invertible with for and for .
and substituting in (1) we obtain
Now suppose that is nonzero for all . We define as
where is the principal logarithm. We also assume that is continuous and differentiable for all . We refer to [Khoa2019, Kliba2019a, Kliba2019d] for the definition of the complex logarithm for a similar change of variables using a high frequency asymptotic behavior for the total field in . Next, this definition was extended in [Khoa2019] to non high values of as long as for those values. To what we know, that asymptotic behavior is not established yet for the two-dimensional case. At the same time, in our numerical studies, we do not see any discontinuity problem with the principal .
We now eliminate by differentiating (8) with respect to
Let be sufficiently large. We approximate the function in (7) and its partial derivative using the truncated Fourier series as
where the coefficients are given by
For each , multiplying both sides of the above equation by and integrating with respect to over , we obtain
Considering two matrices defined as
and an block matrix , each block is an matrix defined as
we can rewrite (12
) as a system of PDEs for the vector valued function
Here the operator is defined as follows: If is an block matrix, each block is an -dimensional column vector and is an -dimensional column vector then is an -dimensional column vector given by
we are able to approximately reformulate the inverse problem as the Cauchy problem for the following system of quasilinear elliptic PDEs:
where and can be computed from the given boundary data and in (3)–(4) using (5), (7) and (11). If we can find by solving problem (14)–(16), the coefficient of interest can be approximately recovered from (8).
We emphasize that the reconstruction algorithm we study in the next section for solving problem (14)–(16) only needs the backscatter data on . In contrast, the convexification method of above cited papers [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019], one has to artificially complete the backscatter data on the other boundaries of for a better stability of computations. On the other hand, the forthcoming analytical results that are mentioned in Introduction for this paper are valid with the Carleman Weight Function (20) only if the Dirichlet data for the system (14) are known on the entire boundary rather than just on its part , i.e. they are valid for those completed data. Thus, our claim in the first sentence of this Remark is based only on our numerical observation and is not supported by the theory. Nevertheless, this numerical observation emphasizes the stability property of our method.
It is well known that the Cauchy problem for an elliptic equation is unstable. Thus, we actually construct a regularization method of solving this problem for our case. A similar numerical method was constructed in [Klib2015] for ill-posed Cauchy problems for a wide class of single quasilinear PDEs, including the elliptic one. However, the Carleman Weight Function used in [Klib2015] for the elliptic case is inconvenient for the numerical implementation since it depends on two large parameters, instead of just one in our case of (20).
3 A numerical reconstruction algorithm
We solve problem (14)–(16) using the weighted quasi-reversibility method. We first make a change of variables to have homogeneous boundary conditions on . Let be a vector valued function which satisfies the boundary conditions (15)–(16). We call the data carrier and its construction is detailed in the numerical study section. Assuming is the solution of problem (14)–(16), we define
Then satisfies the homogeneous boundary conditions on , that is,
Define the function space as
with its associated norm
Let be an arbitrary number. Define the ball as
Next, we define the weighted Tikhonov-like functional as
is a Carleman Weight function and and are constants. From our numerical experience, the regularization terms involving and help us obtain better stability for the computation although they are not needed in the theory of convexification methods in previous studies [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019]. Below we focus on the minimization of the functional on the ball defined in (18). As stated in Introduction, the use of the Carleman Weight Function is inspired by convexification methods whose different versions are described in the above cited publications. For the Carleman Weight Function in (20), a Carleman estimate for the Laplacian has been proved in [Kliba2019c], where the Dirichlet boundary condition is given on the entire boundary of , which requires an artificial complement of the backscatter data given only on the part of . Actually, assigning the Dirichlet data on the entire boundary , one provides an additional stability property to the method. Recall that (see Remark 1) it is our numerical observation that our algorithm only requires the backscatter data on the top boundary of .
An interesting numerical observation is that our algorithm converges and provides better reconstruction results with the Carleman Weight Function, as compared with the case when this function is absent in (19), i.e. when in (20): compare Figure 2 and Figure 3 for numerical results with and without the Carleman Weight Function. The two parameters and along with the regularization parameters , , will be chosen numerically in the implementation of the algorithm. We find the solution as the global minimizer of using a method of gradient descent type. We point out that even though we can prove the global convergence on of the gradient projection method rather than of the gradient descent method, still our numerical observation is that the latter method has good convergence properties. The same observation took place in all previous publications about the convexification where numerical results were presented [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019]. This is a quite useful observation since the numerical implementation of the gradient descent method is much simpler than the one of the gradient projection method.
In this algorithm, we recall that the capital letter notations, for example , are vector valued functions whose components are Fourier coefficients of the corresponding scalar function
with the normal letter notations. We also prescribe a tolerance which forces our iteration to stop after the cost functional no longer decreases much. The tolerance will be chosen numerically in the implementation of the algorithm.
The numerical reconstruction algorithm
Construct the data carrier . Set the initial guess then proceed into the main iteration (Step 2).
Set , compute the cost functional and its gradient .
If and , stop the iteration and move to Step 3.
Otherwise proceed to Step 2b.
Set , where the gradient descent step size is chosen numerically.
Set and compute the corresponding scalar function .
Compute from using the real part of (8) with .
Set and return to Step 2a.
Set , compute from and compute using the real part of (8) with .
We observe from the numerical performance of the algorithm that solving the direct problem to update helps the cost functional decrease faster with respect to iterations. This is important to the algorithm since the cost functional decreases very slowly after the first iteration without this update.
4 Numerical study
In this section, we describe some important details of the numerical implementation of the above algorithm and present some numerical reconstruction results. The first step of the algorithm is to construct the data carrier . Recall that our computational domain . For , define the
and then set
Then for , for , and is a smooth transition from to on . Set
where and are the Cauchy data for , which means on and on . Obviously, these data can be computed from the data and in (3)–(4) using the relation between and in (5) and (7). Then the function satisfies the boundary conditions on and on . Thus, the corresponding vector valued function containing the Fourier coefficients of with respect to that truncated Fourier basis satisfies the boundary conditions (15)–(16). Actually by the definition of the function is zero in , which means that we mainly seek for the scattering object in the upper part of the square since the Cauchy data are given on the top boundary . Indeed, since the Cauchy problem (14)–(16) is unstable, then it is unlikely that even after the regularization, which we do here, one could image scattering objects well, if they are located far from the measurement side . For the parameter we choose in the numerical implementation.
For the implementation of the algorithm, we first discretize the computational domain into uniform grid points , where the mesh size is . The wave number interval is divided into uniform subintervals, where are the midpoints and is the length of each subinterval. Define the lined up index as follows
In this section, using the lined up index we write vector valued functions at grid points as a column vector without changing notations. For instance, for , we have
Let and set and . The weighted Tikhonov-like functional in (19) is discretized using finite differences as
We now describe how to compute for complex valued vector function . Recall that if is a complex variable, is a complex vector and is a complex valued function then we have (see [Kreut2009])
Thus, by the chain rule we have
Treating as a function of complex variables and applying all of the above to each of its summands, we are able to compute and thus obtain using ().
We need a numerical solver for the direct problem (1)–(2) in Step 2b of the reconstruction algorithm and to generate synthetic scattering data for the numerical study. It is well known that the direct problem (1)–(2) is equivalent to the Lippmann-Schwinger integral equation
where is the Hankel function of the first kind of order 0, see [Colto2013]. We exploit the numerical method studied in [Saran2002] to solve this integral equation to generate the Cauchy data and for the inverse problem. Note that the numerical method studied in [Saran2002] assumes smooth coefficients. Its extension to the case of discontinuous coefficients is studied in [Lechl2014] which can be adapted to our discontinuous coefficient examples in this section. We also add an artificial random noise to the data
where is the noise level and are functions taking random complex values and satisfy . We consider noise in the Cauchy backscatter data which means in our numerical examples.
In all numerical examples presented in this section the computational domain is chosen as , where . This means is uniformly discretized into points. The interval of wave numbers is , where . We have found that in (10) is sufficient for the Fourier series truncation, also, see [Khoa2019] for a similar choice. We generate the multifrequency data for the incident plane wave
In the Carleman Weight Function (20), we choose This means that
This choice was made by trial and error, so as choices of all other parameters used in this section. We refer to works on the convexification [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019] for choices of smaller . As to our choice of it seems to give us the optimal results for our numerical examples of this section. The step size of the gradient descent method and the regularization parameters were chosen as:
The tolerance number of the iterations is set up to be . It happens in all our numerical examples that the algorithm stops within 10 iterations. One can observe in the numerical examples that the value of the minimized functional does not change much after 8 or 9 iterations. Since we are interested in , in the final iteration we assign to be zero in the area in which it takes negative values. By our numerical experience, this area is typically below the reconstructed scatterer.
It is important to mention that the initial guess for in all numerical examples below. This goes along well with our theory (to be published) which guarantees that our algorithm converges to the correct solution starting from any point of the ball defined in (18), see item 3 in the end of Introduction. This certainly a significant advantage of our method, compared with locally convergent optimization approaches, which typically need a strong a priori knowledge of the scatterer. Such a knowledge, however, is rarely available in applications.
4.1 Numerical example 1
In this example we consider a single scattering disk characterized by the coefficient which equals 3 inside the disk and zero elsewhere. We can see from the reconstruction result in Figure 2 that the location and and the maximal value of are well reconstructed. It seems to us that the shape of the scattering object is not well-reconstructed because the backscatter data are generated by incident plane waves with a fixed direction, also, see [Nguye2017, Nguye2018, Kliba2019a, Kliba2019d] for similar results. Convergence of the algorithm can be observed from Figure 2(b). From our numerical experience the cost functional does not decrease much after 8 or 9 iterations, see also Figure 2(b).
Now with the numerical result in Figure 3 we want to indicate the importance of the Carleman Weight Function for our numerical algorithm. The algorithm does not converge when the cost functional does not involve the Carleman weight function. Firstly, the error between the cost functionals at two consecutive iterations is never smaller than the tolerance number like what we have when the Carleman Weight Function is present. Therefore, the iterations do not stop with the chosen tolerance. Secondly, the cost functional starts to increase after a certain number of iterations. Figure 3(a) presents the reconstruction result at the sixth iteration where the cost functional obtains its smallest value among 20 iterations. However, this result is not as good as that of Figure 2(c) where the Carleman Weight Function is involved. Indeed, the artifact in Figure 3(a) is slightly stronger and the reconstructed maximal value is 3.2157 while the maximal value of the reconstruction in Figure 2(c) is 3.0014. Also for the next examples, the reconstruction results are always better with the presence of the Carleman Weight Function in the cost functional.
4.2 Numerical example 2
In this example we consider the case of two scattering disks. In the first case in Figure 4(a) two similar scattering disks are considered. The reconstruction result in Figure 4(c) again shows that the algorithm is able to reconstruct very well the location and the maximal values of in this case. The cost functional decreases well within ten iterations, see Figure 4(b). The case of Figure 5(a) is more challenging since the maximal values of in each scattering disk are different. However, the algorithm can provide reasonable reconstruction results in Figures 5(c). One can clearly see the locations of the scattering disks as well as two different maximal values of on each disk. We point out the algorithm in this case can reconstruct the scatterer consisting of two components without using any a priori knowledge about the number of components of the scatterer.
4.3 Numerical example 3
In this example we consider the case of the coefficient which has different values in scattering objects of different shapes. This case is thus more challenging than those of the first two examples. The scatterer in Figure 6(a) consists of a scattering disk in which and a scattering rectangle in which . The reconstruction result in Figure 6(c) shows that the algorithm again can compute the location of the scatterer and the maximal values of in each scattering object. Particularly, we can also see pretty well a difference between the shape of the disk and the rectangle in the reconstruction. The scatterer in Figure 7(a) consists of two scattering disks in which and a scattering rectangle in which . The reconstruction result in Figure 7(c) provides the location and maximal values of the scattering objects. However, the resolution of the reconstruction in this case is not as good as that of the case of two objects since three scattering objects are placed quite close to each other.
We have proposed a new version of the convexification numerical reconstruction method for solving the coefficient inverse scattering problem with multifrequency backscatter data associated with a single direction of the incident plane wave. This method relies on an approximate reformulation of the problem as the Cauchy problem for a system of coupled quasilinear elliptic PDEs. The main ingredients for deriving this formulation are the elimination of the coefficient from Helmholtz equation and the use of truncated Fourier expansion for the total field. To solve the quasilinear elliptic PDE system, we use a weighted quasi-reversibility method in which a Carleman Weight Function is included in the weighted Tikhonov-like functional. The numerical results show that our method is able to efficiently compute the solution without using any a priori information about it. We have shown that values of the dielectric constants of scatterers as well as locations of scatterers can be well reconstructed. However, shapes are not reconstructed well. On the other hand, the recent publication [Khoa2019] shows that all three components of scatterers can be accurately reconstructed in the case when the point source moves along a straight line and frequency is fixed.
Overall, the main advantage of our algorithm, so as other above
cited versions of the convexification method, is its rigorously guaranteed
global convergence, as opposed to the local convergence of the conventional
optimization methods, see Introduction for the definition of the global
convergence. Theoretical analysis of the algorithm as well as its
three-dimensional extension will be addressed in forthcoming publications.
That theoretical analysis will consists in detailed proofs of results
announced in items 1-3 of Introduction.
Acknowledgement. The authors are grateful to Loc Nguyen and Khoa Vo from the University of North Carolina at Charlotte for helpful discussions. The first and second authors are partially supported by NSF grant DMS-1812693. The third author is supported by US Army Research Laboratory and US Army Research Office grant W911NF-19-1-0044.