 # A numerical reconstruction algorithm for the inverse scattering problem with backscatter data

This paper is concerned with the inverse scattering problem which aims to determine the coefficient of the Helmholtz equation from multifrequency backscatter data. We propose an efficient numerical algorithm to solve this nonlinear and ill-posed inverse problem without using any advanced a priori knowledge of the solution. To study the algorithm we first eliminate the coefficient from the Helmholtz equation using a change of variables. Then using a truncated Fourier expansion for the wave field we approximately reformulate the inverse problem as a system of quasilinear elliptic PDEs, which can be numerically solved by a quasi-reversibility approach. The cost functional for the quasi-reversibility method is constructed as a Tikhonov-like functional that involves a Carleman weight function. Our numerical study shows that using a method of gradient descent type one can find the minimizer of this Tikhonov-like functional without any advanced a priori knowledge about it.

## Authors

##### 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

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 function with 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 [Kliba2019a, Kliba2017, Kliba2019c, Khoa2019] for a similar assumption. Suppose that the object is illuminated by some incident plane wave . 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

 Δu+k2(1+a(x))u=0,x∈R2, (1) lim|x|→∞|x|(∂(u−uin)∂|x|−ik(u−uin))=0. (2)

Here is the wavenumber. 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

 Ω=(−R,R)2,Γ=(−R,R)×{R}. (3)

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

 g0(x,k) :=u(x,k),for x∈Γ,k∈[k––,¯¯¯k], (4) g1(x,k) :=∂u∂x2(x,k),for x∈Γ,k∈[k––,¯¯¯k], (5)

where is the total wave in the scattering problem (1)–(2). Determine the material parameter in (1) for , see also Figure 1 for a schematic of the inverse problem.

This inverse problem belongs to the wider class of inverse scattering problems which in general aim to recover information about the coefficient (e.g. its domain of definition and/or its values) from 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. We only discuss some of the numerical solution methods for the interest of the present paper. The most studied approach is probably optimization based methods, see for instance

[Bakus2004, Engl1996]. However, it is well known that these methods may suffer from local minima and their convergence analysis is also not known 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, only geometrical information of the scatterer can be reconstructed with qualitative methods. Furthermore, these methods typically require muti-static data which are sometime not available in some practical applications. Figure 1: Schematic of the inverse scattering from a penetrable bounded object characterized by the function a(x). The incident plane wave propagates downward toward the scattering object. The backscatter data are measured on the top boundary Γ of the computational domain Ω.

The method studied in this paper is an extended study from a recently new approach called globally convergent numerical methods (GCNM) for solving coefficient inverse problems. The GCNM typically aims to reconstruct a coefficient in inverse scattering problems using scattering data for a single direction of the incident plane wave. Its main advantage is to avoid the local minimum problem suffered by optimization based methods. 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. Our method is inspired by the second type of the GCNM that has been more recently studied in [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019]. The central idea of the GCNM of the second type is the construction of globally convex functionals involving Carleman weight functions. The first step is to eliminate the coefficient from the Helmholtz equation using a change of variables. Using a truncated Fourier expansion for the total field we approximately reformulate the inverse problem as a system of quasilinear elliptic PDEs. We then propose a quasi-reversibility method to solve the elliptic PDE system. Inspired by the GCNM of the second type mentioned above the cost functional in the quasi-reversibility method is involved with a Carleman weight function which plays an important role in the numerical performance of the method. A method of gradient descent type is exploited to find the global minimizer of the cost functional without using any advanced a priori information about it.

Comparing with the recent works of the GCNM of the second type cited above, 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. This work is mostly related to [Kliba2018a] in which the one-dimensional version of the inverse problem has been studied. Using the new Fourier basis from [Kliba2017] the multi-dimensional 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 iteration. Thirdly, unlike the previous works in [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 constant functions. 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 will be addressed in an incoming publication.

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 formulate our inverse problem as 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]

 ψn(k)=(k−k0)n−1ek−k0,k∈(k––,¯¯¯k),n=1,2,…

Applying the Gram–Schmidt process to

we obtain an orthonormal basis which has the following properties, see also [Kliba2017]:
i) for all
ii) The matrix , where and , is invertible with for and for .

Now setting

 p(x,k)=u(x,k)uin(x,k), (6)

and substituting in (1) we obtain

 Δp(x,k)+k2a(x)p(x,k)−2ik∂x2p(x,k)=0. (7)

Now suppose that is nonzero for all . We define as

 v(x,k)=log(p(x,k))k2, (8)

where is the principal logarithm. We also assume that is continuous and differentiable for all . We refer to [Khoa2019] for the definition of the complex logarithm for a similar change of variables using a high frequency asymptotic behavior for the total field in . To what we know this asymptotic behavior is not established yet for the two-dimensional case. We do not see any discontinuity problem with the principal in the numerical implementation. Using (8) we substitute in (7) and rewrite (7) in terms of as follows

 Δv(x,k)+k2∇v(x,k)⋅∇v(x,k)−2ik∂x2v(x,k)+a(x)=0. (9)

We now eliminate by differentiating (9) with respect to

 Δ(∂kv)+2k∇v⋅∇(v+k∂kv)−2i(∂x2v+k∂x2∂kv)=0. (10)

Now let be sufficiently large. We approximate the function in (8) and its partial derivative using the truncated Fourier series as

 v(x,k)=N∑n=1vn(x)Φn(k),∂kv(x,k)=N∑n=1vn(x)Φ′n(k)

where the coefficients are given by

 vn(x)=∫¯¯¯kk––v(x,k)Φn(k)dk. (11)

Using these truncated series we approximate (10) by

 N∑n=1Φ′n(k)Δvn(x)+2kN∑n=1N∑l=1Φn(k)(Φl(k)+kΦ′l(k))∇vn(x)⋅∇vl(x) −2iN∑n=1(Φn(k)+kΦ′n(k))∂x2vn(x)=0.

For each , multiplying both sides of the above equation by and integrating with respect to over , we obtain

 N∑n=1(∫¯¯¯kk––Φm(k)Φ′n(k)dk)Δvn(x)+N∑n=1N∑l=1(2k∫¯¯¯kk––Φm(k)Φn(k)[Φl(k)+kΦ′l(k)]dk)∇vn(x)⋅∇vl(x)−N∑n=1(2i∫¯¯¯kk––Φm(k)[Φn(k)+kΦ′n(k)]dk)∂x2vn(x)=0. (12)

Considering two matrices defined as

 D =(dmn),dmn=∫¯¯¯kk––Φm(k)Φ′n(k)dk, S =(smn),smn=−2i∫¯¯¯kk––Φm(k)[Φn(k)+kΦ′n(k)]dk,

and an block matrix , each block is an matrix defined as

 b(l)mn=2k∫¯¯¯kk––Φm(k)Φn(k)[Φl(k)+kΦ′l(k)]dk,

we can rewrite (12

) as a system of PDEs for the vector valued function

 DΔV(x)+B∂x1V(x)∙∂x1V(x)+B∂x2V(x)∙∂x2V(x)+S∂x2V(x)=0. (13)

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

 P∙V=⎡⎢ ⎢ ⎢ ⎢⎣P1⋅VP2⋅V⋮PN⋅V⎤⎥ ⎥ ⎥ ⎥⎦.

Defining

 Q(V)=DΔV+B∂x1V∙∂x1V+B∂x2V∙∂x2V+S∂x2V,

we are able to approximately reformulate the inverse problem as the following system of quasilinear elliptic PDEs

 Q(V) =0in Ω, (14) V =G0on Γ, (15) ∂x2V =G1on Γ, (16)

where and can be computed from the given boundary data and in (4)–(5) using (6),(8) and (11). If we can find by solving problem (14)–(16), the coefficient of interest can be approximately recovered from (9).

###### Remark 1

We want to emphasize that the reconstruction algorithm we study in the next section for solving problem (14)–(16) only needs the backscatter data on . In the convexification method papers [Kliba2019a, Kliba2017, Kliba2019c, Khoa2019] one has to artificially complete the backscatter data on the other boundaries of for a better stability of computations.

## 3 A numerical reconstruction algorithm

We solve problem (14)–(16) using 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

 W=V−F.

Then satisfies the homogeneous boundary conditions on , that is,

 W=∂x2W=0on Γ.

We aim to find in the following function space

 X={W∈[H2(Ω)]N, W=∂x2W=0 on Γ}

with its associated norm

 ∥W∥X=(N∑n=1∥wn∥2H2(Ω))12,where W(x)=[w1(x) w2(x) … wN(x)]T.

Next we define the cost functional

 J(W)=∫Ω|Q(W+F)|2φ2dx+ρ∥W∥2X+α1∫Γ|W|2dx+α2∫Γ∣∣∂x2W∣∣2dx, (17)

where

 φ(x)=e−λ(x2−s)2 (18)

is a Carleman weight function and and are constants. The use of the Carleman weight function is inspired by convexification methods recently developed in [Kliba2018a, Kliba2019a, Kliba2017, Kliba2019c, Kliba2019d, Khoa2019]. In these papers, the authors studied globally convex functionals for solving nonlinear coefficient inverse problems with the help of Carleman weight functions. For the Carleman weight function in (18

), a Carleman estimate for the Laplacian has been proved in

[Kliba2018a] where the scattering data are given on the entire boundary of . Recall that (see Remark 1) our algorithm only requires the backscatter data on the top boundary of .

We observe from the numerical simulations that the algorithm converges and provides better reconstruction results with the Carleman weight function, see Figure 2 and Figure 3 for numerical results with and without the Carleman weight function. The two constants 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. In the following we will describe our numerical algorithm for finding the coefficient for the inverse problem in (4)–(5), in which finding the minimizer of is one of the main components of the algorithm.

###### Remark 2

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 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

#### Step 1.

Construct the data carrier . Set the initial guess then proceed into the main iteration (Step 2).

#### Step 2a.

Set , compute the cost functional and its gradient .

• If and , stop the iteration and move to Step 3.

• Otherwise proceed to Step 2b.

#### Step 2b.

• Set , where the gradient descent step size is chosen numerically.

• Set and compute the corresponding scalar function .

• Compute from using (9) with .

• Compute by solving the direct problem (1)–(2) with , and set .

• Compute from using (6) and (8), and then compute .

#### Step 3.

Set , compute from and compute using (9) with .

###### Remark 3

We observe from the numerical performance of the algorithm that solving the direct problem to update helps the cost functional decrease faster in the iteration. 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 about the numerical implementation of the algorithm and present some numerical reconstruction results obtained from the algorithm. The first step of the algorithm is to construct the data carrier . Recall that our computational domain . For , define

 χ0(t)={exp(−Rt+ξ),t>−ξ0,t≤−ξ

and

 χ(t)=χ0(t)χ0(t)+χ0(R−t−2ξ).

Then for , for , and is a smooth transition from to on . Define

 f(x,k)=[~g0(x,k)+(x2−R)~g1(x,k)]χ(x2),

where and are computed from the data and using (6) and (8). Then satisfies the boundary conditions on and on . Thus the corresponding vector valued function containing the Fourier coefficients of with respect to the 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 in which the Cauchy data are given on the top boundary . 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

 m=m(i,j,r)=i+(j−1)(Nx+1)+(r−1)(Nx+1)2,1≤i,j≤Nx+1, 1≤r≤N.

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

 U=[um],1≤m≤(Nx+1)2N,

where

 um=um(i,j,r)=ur(xij).

Let and set and . The cost functional in (17) is discretized using finite differences as

 J(W) +N∑s=1bsmrh2x(^wm(i,j+1,r)−^wm(i,j,r))(^wm(i,j+1,s)−^wm(i,j,s)) +N∑s=1bsmrh2x(^wm(i+1,j,r)−^wm(i,j,r))(^wm(i+1,j,s)−^wm(i,j,s)) +smrhx(^wm(i+1,j,r)−^wm(i,j,r))]φ(xij)∣∣∣2 +ρh2xN∑m=1Nx+1∑j=1Nx+1∑i=1∣∣wm(i,j,m)∣∣2 +ρh2xN∑m=1Nx∑j=2Nx∑i=2[∣∣∣wm(i,j+1,m)−wm(i,j,m)hx∣∣∣2+∣∣∣wm(i+1,j,m)−wm(i,j,m)hx∣∣∣2 +∣∣∣wm(i,j+1,m)−2wm(i,j,m)+wm(i,j−1,m)h2x∣∣∣2+∣∣∣wm(i+1,j,m)−2wm(i,j,m)+wm(i−1,j,m)h2x∣∣∣2 +α1hxN∑m=1Nx+1∑j=1∣∣wm(Nx+1,j,m)∣∣2+α2hxN∑m=1Nx∑j=2∣∣∣wm(Nx+1,j,m)−wm(Nx,j,m)hx∣∣∣2.

Now we describe how to compute for complex variable . Recall that if is a complex variable, is a complex vector and is a complex-valued function then we have (see [Kreut2009])

 i.∂∂z|z|2=∂∂z(z¯z)=¯z, ii.∂h∂z(z)=[∂h∂z1 ∂h∂z2 … ∂h∂zM], iii.∇h(z)=(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂h∂z(z))T.

Thus by the chain rule we have

 ∂∂z|h|2(z)=¯¯¯¯¯¯¯¯¯¯h(z) ∂h∂z(z).

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

 u(x,k)=uin(x,k)+k2∫Ωi4H(1)0(k|x−y|)a(y)u(y,k)dy, (19)

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 artificial noise to the these data

 gj(x,k) =gj(x,k)+δ∥gj∥L2Nj(x,k),j=0,1,

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.

Now for the numerical examples presented in this section the computational domain is chosen as , where . This means is uniformly discretized into points. The wave number interval , where . We find that is sufficient for the Fourier series truncation, see also [Khoa2019] for a similar choice. We generate the multifrequency data for the incident plane wave

 uin(x)=e−ikx2,k∈[0.5,2].

For the Carleman weight function, we choose which means

 φ(x)=e−5(x2−1)2.

This choice is made by test and error. We refer to [Khoa2019, Kliba2019a] and related works on convexification methods for choices of smaller . This choice of seems to give us the optimal results for our numerical examples in this section. The the gradient descent step size and the regularization parameters are also determined by test and error as

 ε=α1=10−3,ρ=α2=10−5.

and the tolerance in the iteration is set up to be . It happens in all the examples that the algorithm stops within 10 iterations. As one can see in the numerical examples that the cost 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 it takes negative values. From our numerical experience this area is typically below the reconstructed scatterer.

### 4.1 Numerical example 1

In this example we consider a single scattering ball characterized by the coefficient which equals 3 inside the ball 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, see also [Nguye2017, Nguye2018, Kliba2019a] 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 seem to 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 like what we have with the presence of the Carleman weight function. Therefore, the iterations does not stop with the chosen tolerance. Secondly, the cost functional starts to increase after a 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 there is the Carleman weight function. 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 balls. In the first case in Figure 4(a) two similar scattering balls 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 ball are different. However, the algorithm can provide reasonable reconstruction results in Figures 5(c). One can clearly see the location of the scattering balls as well as two different maximal values of on each ball. We want to mention that 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. We can see that the initial guess for is just zero in all numerical examples. This can be an advantage comparing with local optimization approaches which typically need strong a priori knowledge of the scatterer.

### 4.3 Numerical example 3

In this example we consider the case of the coefficient which have 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 ball 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 ball and the rectangle in the reconstruction. The scatterer in Figure 7(a) consists of two scattering balls 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.