Convexification method for nonlinear SAR imagingMichael V. Klibanov et al.

# Convexification Inversion Method for Nonlinear SAR Imaging with
Experimentally Collected Data ††thanks: Submitted to the editors DATE. The
work of Klibanov, Khoa, Smirnov, Nguyen, Bidney and Astratov was supported by
US Army Research Laboratory and US Army Research Office grant
W911NF-19-1-0044.

Michael V. Klibanov Corresponding author. Department of Mathematics and
Statistics, University of North Carolina at Charlotte, Charlotte, NC, 28223
USA<EMAIL_ADDRESS>Vo Anh Khoa Department of Mathematics and
Statistics, University of North Carolina at Charlotte, Charlotte, NC, 28223
USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Alexey V. Smirnov Department
of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L
3G1<EMAIL_ADDRESS>Loc H. Nguyen33footnotemark: 3 Grant W. Bidney
Department of Physics and Optical Science, University of North Carolina at
Charlotte, Charlotte, NC 28223, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Lam H. Nguyen U.S. Army Research Laboratory, Adelphi, MD 20783-1197, USA (, ).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Anders J.
Sullivan66footnotemark: 6 Vasily N. Astratov55footnotemark: 5

###### Abstract

This paper is concerned with the study of a version of the globally convergent
convexification method with direct application to synthetic aperture radar
(SAR) imaging. Results of numerial testing are presented for experimentally
collected data for a fake landmine. The SAR imaging technique is a common tool
used to create maps of parts of the surface of the Earth or other planets.
Recently, it has been applied in the context of non-invasive inspections of
buildings in military and civilian services. Nowadays, any SAR imaging
software is based on the Born approximation, which is a linearization of the
original wave-like partial differential equation. One of the essential
assumptions this linearization procedure needs is that only those dielectric
constants are imaged whose values are close to the constant background. In
this work, we propose a radically new idea: to work without any linearization
while still using the same data as the conventional SAR imaging technique
uses. We construct a 2D image of the dielectric constant function using a
number of 1D images of this function obtained via solving a 1D coefficient
inverse problem (CIP) for a hyperbolic equation. Different from our previous
studies on the convexification method with concentration on the global
convergence of the gradient projection method, this time we prove the global
convergence of the gradient descent method, which is easier to implement
numerically.

###### keywords:

SAR imaging, gradient descent method, convexification, global convergence,
experimental data, data propagation

78A46, 65L70, 65C20

##  1 Introduction

SAR imaging is a popular technique broadly used by aircraft, drones and
satellites to map some parts of the surface of planets including Earth; cf.
e.g. [4, 8, 12, 33]. The SAR imaging technique exists for several years and is
developed so well by now that those flying devices commonly use a commercial
SAR software. The Physics behind it can be understood in the following manner.
It is known that the size of an antenna of any flying device is small.
However, if combining those antennas running along a large distance covered by
that device, then the size of that synthetic antenna becomes quite large. As a
result, one obtains a good resolution of SAR images.

According to [12] and references cited therein, the foundational assumption of
any SAR imaging software is the Born approximation. From the mathematical
standpoint, this means that the dielectric constant of the medium is assumed
to be a small perturbation of the unity. In other words, the corresponding
Coefficient Inverse Problem (CIP) for a wave-like hyperbolic partial
differential equation (PDE) is linearized near that unity. This assumption
essentially enables one to obtain an explicit formula for the solution of that
CIP, which is undoubtedly attractive from numerical standpoint. Up-to-date,
any commercial software for SAR imaging relies on that formula as well as some
of its variants in different circumstances of the above-mentioned
applications. In our first work of this direction [21], we have demonstrated
numerically that the Born approximation leads to significant errors in values
of the dielectric constants of targets.

It is worth mentioning that the motivation of [21] about nonlinear SAR imaging
stems from the non-invasive inspections of buildings with a particular
concentration on the identification of improvised explosive devices (IEDs);
cf. [2, 4, 26]. Needless to say, having knowledge of the dielectric constant
of the medium in that perspective is an important piece of information which
hopefully might help to decrease the false alarm rate for military purposes.
With this application being in mind, there have been numerous models with
distinctive numerical approaches tackled by the first author and his research
team for almost a decade; cf. e.g. [6, 14, 16, 18, 22, 23, 25, 39]. The
interests therein dwell not only in the approximations of the dielectric
constants, but in locations, sizes and even shapes of such targets. This
information is anticipated to be helpful in the development of future
classification algorithms that would better distinguish between explosives and
clutter in the battlefield. While the conventional SAR imaging might only be
useful in detecting the object’s location and shape, our wish here is to fill
the gap regarding the identification problem of the material of the
constituent medium. Except of our work [21], we are unaware of any
publication, where the dielectric constants would be indeed computed for SAR
data, albeit the theory of SAR imaging tells one that those should be imaged
[12].

The work [16] is our first attempt to develop our convexification method (see
this section below about this method) for a CIP for the Helmholtz equation
with the moving point source. Next, we have shown that the technique of [16]
works well on experimentally collected data [14, 15]. The case of the moving
point source is somewhat similar to the case of SAR imaging. We however remark
that the SAR technique itself follows time-dependent models rather than what
we have done so far in [14, 15, 16] with the frequency domain case with a
fixed frequency. Yet, another difference between those recent works and the
SAR scenario is that in the SAR case the incident waves are sent from a
specific altitude angle and are radiated not by a point source but originated
in a 2-D area, which acts as an antenna. Therefore, a _radically new_ idea has
to be proposed to overcome these difficulties.

We also note that various CIPs with the moving source and a fixed frequency
were considered by the group of R. G. Novikov since 1988 [28], see [1, 3, 27,
29]. However, statements of CIPs in those publications are different from ones
in [14, 15, 16] and reconstruction procedures are also different from the
convexification.

The CIP of our interest with the SAR data is actually an underdetermined one
[21]. Indeed, the sought function that characterizes the spatial distribution
of the dielectric constant depends on three variables, while the SAR data
depemd on two variables. It is the underdertermination of the SAR data which
causes the _main challenge_ of the entire effort. To address this, we form a
2D image via solutions of a number of 1D CIPs for a wave-like hyperbolic
equation. Following the approach commenced in [37], we have established in
[21] a transformation of the underlying hyperbolic equation into a nonlocal
nonlinear PDE. Then the solution is approximated using the minimization of a
weighted globally strictly convex Tikhonov-like functional.

As to the minimization problem, we emphasize that we solve that 1D CIP by a
version of the convexification method, which is presented in [37]. The
convexification is a concept that the research team of the first author has
been pursuing for a number of years. This concept is based on the
Bukhgeim–Klibanov method [7], which was originally introduced in 1981 only for
proofs of global uniqueness theorems for CIPs and without any idea at that
time of its significant potential in applications to numerical methods for
CIPs. The key ingredient of the convexification method is the construction of
a Tikhonov-like cost functional weighted by a suitable Carleman weight
function. The presence of the Carleman weight function is needed to
“convexify” the cost functional involving the corresponding PDE operator.
Driven by the so-called Carleman estimate, this method is well known to be a
_globally convergent_ numerical method for CIPs. We call a numerical method
for a CIP globally convergent if there is a theorem which claims that this
method delivers at least one point in a sufficiently small neighborhood of the
true solution without any advanced knowledge of this neighborhood. In the past
works, the convexification method and its variants were designed due to the
need of solving CIPs that involves nonlinear differential operators. It is
obvious that for nonlinear cases, the conventional least squares method can be
trapped in local minima and ravines; cf. e.g. [34, Figure 1] for a typical
graphical observation when local minima occur. In the convexification, we are
capable to prove the strict convexity of the functional on a given convex
bounded set of a Hilbert space. The global convergence is done when one can
prove the convergence of the scheme towards the true solution starting from
any element in that bounded set, whose diameter is arbitrary.

As any CIP is highly nonlinear, we usually call our technique “nonlinear SAR
imaging”. We point out that our method is a heuristic one since we cannot
derive it from Physics. We can justify our idea only numerically. At the same
time, we note that the conventional SAR imaging also uses a number of
heuristic approximations (in which the Born approximation is the first one) to
finalize their explicit formula. The only rigorous part of our method is the
global convergence theory of the convexification method for our CIP.

Although the essence of this theory is taken from our recent work [37], we
still have a significantly new analytical element here. More precisely, in all
previous works on the convexification (e.g., [14, 16, 19, 20, 22, 36, 37]) the
global convergence was established only for the gradient projection method,
which is much harder to implement numerically than the regular gradient
descent method. On the other hand, here we establish, _for the first time_ ,
that actually the gradient descent method being applied to our weighted
Tikhonov-like cost functional converges globally. This result explains our
consistently repeated observation of all previous studies of convexification,
in which only the gradient descent method was numerically implemented and
numerical results were good ones; see, e.g. [5, 14, 15, 16, 18, 20, 21, 36,
37]. Hence, the extended theory in this work is also a new mathematical aspect
of our first publication [21] which is a more engineering one. Compared with
[21], in this paper we present new numerical results for nonlinear SAR
imaging. Numerical results are presented for experimentally collected SAR
data. Experimental data were collected at the University of North Carolina at
Charlotte.

The paper is organized as follows. The contemporary setting of SAR imaging is
shortly introduced in section 2. Then section 3 is devoted to the statements
of the forward and inverse problems. We also provide an asymptotic analysis of
the incident wave of the forward problem in Lemma 1. This is also a novelty of
the previous work [21] from the theoretical standpoint. Following [37], we
revisit the formulation of the 1D CIP in section 4. Here, we also present the
convexification method for this CIP. Our main results are then in section 4.3
for the global convergence of the gradient descent method; cf. Theorems 7 and
8. As soon as all mathematical apparatus is completed, we work in section 5
with our experimental results. In that place, we prove experimentally that our
proposed method delivers accurate values of the dielectric constant of an
explosive-like target. In this regard, the experimental data are those
backscattering ones for the experimental target buried in the sandbox.
Finally, some conclusions follow in section 6.

##  2 Settings of SAR imaging

We consider in this paper the simplest version of the SAR imaging, called
stripmap SAR geometry (cf. e.g. [35]). In the context of the inspection of
buildings, the SAR device consists of an antenna and a receiver. They move
concurrently along a straight line and cover a certain finite interval of this
line. Due to the large antenna/targets distances, we assume below, when
discussing the inverse problem only, that the radiating antenna and the
receiver coincide and actually, this is a point. This is a typical assumption
of SAR imaging; cf. [12]. While moving, the antenna radiates pulses of a time
resolved component of the electric wave field, and the receiver collects the
backscattering time resolved signal. Even though the propagation of the
electromagnetic wave field is governed by Maxwell’s equations, the SAR imaging
works with only a single wave-like PDE governing the propagation of a single
component of the electric wave field [12].

Let $\Omega=(-\overline{R},\overline{R})^{3}\subset\mathbb{R}^{3}$ for
$\overline{R}>0$ be a cube in which our targets of interest dwell. Given
$y_{0},z_{0}\in\mathbb{R}$, we denote
$\mathbf{x}_{0}=(x_{0},y_{0},z_{0})\in\mathbb{R}^{3}$ by the single point,
where the transmitter and the receiver are assumed to coincide. Henceforth, by
letting $x_{0}$ run along an interval $(-L,L)$ for $L>0$, we define the line
of sources as follows:

$L_{\text{src}}:=\left\\{\mathbf{x}_{0}=\left(x_{0},y_{0},z_{0}\right)\in\mathbb{R}^{3}:x_{0}\in\left(-L,L\right)\right\\},$

where numbers $y_{0}$ and $z_{0}$ are fixed. Consider the straight line
passing through the point $\mathbf{x}_{0},$ orthogonal to the line of sources
$L_{\text{src}}$ and forming the angle $\theta\in\left(0,\pi/2\right)$ with
the plane $\left\\{z=z_{0}\right\\}.$ The half $CRL(\mathbf{x}_{0})$ of this
straight line with the end point $\mathbf{x}_{0},$ which looks towards the
domain of our interest $\Omega$, is called “the central line” of the antenna.
The angle $\theta$ remains the same for all positions of the source, is called
the “elevation angle” and it characterizes the direction of the incident wave
field. Thus, central lines of all antennas are parallel to each other.
Consider the plane $P$ which is passing through all these central lines as
well as through the line of sources. $P$ is called “slant range” or “imaging
plane”. We remark that in the context of SAR imaging, the domain of interest
$\Omega$ should be practically located far from the line of source
$L_{\text{src}}$.

These settings of the geometry of the SAR imaging are depicted in Figure 1. We
have prepared a glass bottle filled with clear water for our experiment; see
Figure 6a for some details about this item. As to the application being in
mind, this bottle is well-suited for mimicking the usual Glassmine 43 ([30]),
a non-metallic anti-personnel land mine having a glass body. It was used in
the World War II era when the Germans wanted to obstruct the detection service
of the Allies; cf. [30].

In order to work with the inverse problem of SAR imaging, we pose the forward
problem in the next section. This step helps us to understand the qualitative
behavior of the wave field, which is incident by a close to realistic antenna,
see Lemma 1 and our comment below its proof.

Figure 1: A schematic diagram for a stripmap SAR imaging. This depiction
follows our experimental configuration in section 5. Buried in a wooden-like
box filled with the dry sand, the experimental target is actually a bottle of
water mimicking a fake landmine. Here the transmitter and receiver run along a
straight line in the $x-$direction. Their positions coincide and are denoted
as sources $\mathbf{x}_{0}$ (yellow dots). Our idea is that the 2D image on
the slant range (red dashed imaging plane) is formed by combining parallel 1D
images on that plane. Those one-dimensional images are formed via solutions to
the CIP stated in section 4.1.

##  3 Statements of forward and inverse problems

###  3.1 Forward problem

For $\mathbf{x}=\left(x,y,z\right)\in\mathbb{R}^{3}$, consider a smooth
function $\varepsilon_{r}(\mathbf{x})$ that represents the spatially
distributed dielectric constant of the medium. We assume that

(1) $\begin{cases}\varepsilon_{r}\left(\mathbf{x}\right)\geq 1&\text{for
}\mathbf{x}\in\Omega,\\\ \varepsilon_{r}\left(\mathbf{x}\right)=1&\text{for
}\mathbf{x}\in\mathbb{R}^{3}\backslash\Omega.\end{cases}$

In (1), $\varepsilon_{r}\left(\mathbf{x}\right)=1$ means that we have vacuum
outside of the domain of interest.

Now, let $T>0$ and $\mathbf{x}_{0}\in L_{\text{src}}$. Cf. [12], we work below
with the following initial value (or Cauchy) problem:

(2) $\begin{cases}\varepsilon_{r}\left(\mathbf{x}\right)u_{tt}=\Delta
u+W\left(t,\mathbf{x},\mathbf{x}_{0}\right)&\text{for
}\mathbf{x}\in\mathbb{R}^{3},t\in\left(0,T\right),\\\
u\left(\mathbf{x},\mathbf{x}_{0},0\right)=u_{t}\left(\mathbf{x},\mathbf{x}_{0},0\right)=0&\text{for
}\mathbf{x}\in\mathbb{R}^{3}.\end{cases}$

Here, $u=u(\mathbf{x},\mathbf{x}_{0},t)$ is the amplitude of the time resolved
component of the electric field, and the source function $W$ is defined as:

(3)
$W\left(t,\theta,\mathbf{x}-\mathbf{x_{0}}\right)=A\left(t\right)e^{-\text{i}\omega_{0}t}\widetilde{W}\left(\theta,\mathbf{x}-\mathbf{x}_{0}\right),\quad
A(t)=\chi_{\tau}(t)e^{-\text{i}\alpha\left(t-\tau/2\right)^{2}},$

where $\omega_{0}$ is the carrier frequency, $\alpha$ is the chirp rate, and
$\tau$ is its duration. In SAR imaging, the expression (3) is called “linear
modulated pulse” or “chirp” with

$\chi_{\tau}\left(t\right)=\begin{cases}1&\text{for
}t\in\left(0,\tau\right),\\\ 0&\text{for }t<0.\end{cases}$

Different from what we have assumed above that the antenna is a point source,
when working with the data generation, we take into account a more realistic
scenario: to use a regular circular antenna. This is presented in our
establishment of the function $\widetilde{W}$ in (3) as follows: The circular
antenna is a disk of the diameter $D>0$ in the plane $P_{\text{ant}}$. Denoted
by $S(\mathbf{x}_{0},\theta,D)$, this disk is centered at $\mathbf{x}_{0}$.
Now, $\mathbf{x}_{0}$ lies in the line $L_{\text{scr}}$ that intersects with
the plane $P_{\text{ant}}$. On the other hand, $P_{\text{ant}}$ is orthogonal
to the central line $CRL(\mathbf{x}_{0})$. Next, we change coordinates via a
rotation
$\mathbf{x}\rightarrow\mathbf{x}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$
such that the central line $CRL(\mathbf{x}_{0})$ of the antenna is parallel to
the $z^{\prime}$-axis. Then
$\mathbf{x}_{0}=(x_{0},y_{0}^{\prime},z_{0}^{\prime})\in L_{\text{scr}}$.

Given a number $\eta\in(0,D/2)$, we introduce a smooth function
$m(\eta,\theta,\mathbf{x},\mathbf{x}_{0},D)$:

(4) $\displaystyle
m\left(\eta,\theta,\mathbf{x},\mathbf{x}_{0},D\right)=\begin{cases}1&\text{for
}\left|\mathbf{x}-\mathbf{x}_{0}\right|<D/2-\eta,\\\ 0&\text{for
}\left|\mathbf{x}-\mathbf{x}_{0}\right|\geq D/2,\\\
\in\left[0,1\right]&\text{for
}D/2-\eta\leq\left|\mathbf{x}-\mathbf{x}_{0}\right|<D/2.\end{cases}$

In this regard, $m(\eta,\theta,\mathbf{x},\mathbf{x}_{0},D)=0$ outside of the
ball centered at $\mathbf{x}_{0}$ with the radius $D/2$. Henceforth, we define
the function $\widetilde{W}$ in (3) as

(5)
$\displaystyle\widetilde{W}(\theta,\mathbf{x}-\mathbf{x}_{0})=\delta(z^{\prime}-z_{0}^{\prime})m(\eta,\theta,\mathbf{x}-\mathbf{x}_{0},D),$

where $\delta$ is the so-called delta function. By this way, due to the
presence of $m$ in (4), the source function $W$ vanishes outside of our disk
antenna $S(\mathbf{x}_{0},\theta,D)$.

######  Remark 1.

We assume that the linear size of the antenna is small compared with the
distance between that antenna and the domain of interest $\Omega$. This is
equivalent with the assumption that the target is located far from the
antenna. In physical optics, the far-field region is conceptually determined
by the Fresnel criterion. Cf. [12, Chapter 2], the full angular beamwidth of
the underlying antenna is $\Psi=2\lambda_{0}/D$, where $\lambda_{0}=(2\pi
c_{0})/\omega_{0}$ is the wavelength of the signal of the antenna and $c_{0}$
is the speed of light in the vacuum. Such information is only used in the SAR
data generation as our best attempt to model the realistic scenario.

The SAR data is actually the function
$F(\mathbf{x}_{0},t)=u(\mathbf{x}_{0},\mathbf{x}_{0},t)$ for
$\mathbf{x}_{0}\in L_{\text{scr}}$, $t>0$. We generate these data via solving
the forward problem (2). When doing so, we apply the Fourier transform with
respect to the time variable $t$, which reads as

(6)
$v(\mathbf{x},\mathbf{x}_{0},k)=\int_{0}^{\infty}u(\mathbf{x},\mathbf{x}_{0},t)e^{ikt}dt,$

assuming that the integral converges. Then, it is straightforward to obtain
the Helmholtz equation for the resulting function $v$. Subsequently, in view
of the fact that $\varepsilon_{r}(\mathbf{x})=1$ in
$\mathbb{R}^{3}\backslash\Omega$, we are led to an analog of the so-called
Lippmann–Schwinger equation (cf. e.g. [9]):

(7)
$v\left(\mathbf{x},\mathbf{x}_{0},k\right)=v_{0}\left(\mathbf{x},\mathbf{x}_{0},k\right)+k^{2}\mathop{\displaystyle\int}\limits_{\Omega}\frac{\text{exp}\left(\text{i}k\left|\mathbf{x}-\mathbf{x}^{\prime}\right|\right)}{4\pi\left|\mathbf{x}-\mathbf{x}^{\prime}\right|}\left(\varepsilon_{r}\left(\mathbf{x}^{\prime}\right)-1\right)v\left(\mathbf{x}^{\prime},\mathbf{x}_{0},k\right)d\mathbf{x}^{\prime}.$

Here, the incident wave $v_{0}\left(\mathbf{x},\mathbf{x}_{0},k\right)$ has
the following form:

$\displaystyle
v_{0}\left(\mathbf{x},\mathbf{x}_{0},k\right):=\mathop{\displaystyle\int}\limits_{0}^{\tau}A\left(t\right)e^{-\text{i}\omega_{0}t}e^{\text{i}kt}dt$
$\displaystyle\times\mathop{\displaystyle\int}\limits_{S\left(\mathbf{x}_{0},\theta,D\right)}\frac{\text{exp}\left(\text{i}k\sqrt{\left(x-\xi_{1}\right)^{2}+\left(y-\xi_{2}\right)^{2}+\left(z-z_{0}\right)^{2}}\right)}{4\pi\sqrt{\left(x-\xi_{1}\right)^{2}+\left(y-\xi_{2}\right)^{2}+\left(z-z_{0}\right)^{2}}}m\left(\eta,\theta,\boldsymbol{\xi},\mathbf{x}_{0},D\right)d\xi_{1}d\xi_{2}.$

In (7), $\text{i}=\sqrt{-1}$. The number $k>0$ is known as the wavenumber and
we will solve (7) for $k\in[k_{1},k_{2}]$, where $k_{1}$ and $k_{2}$ are
chosen numerically in our simulations. Thus, applying the inverse Fourier
transform to the function $v$ with the integration over that interval
$[k_{1},k_{2}]$ to (6) would result in our targeted SAR data
$F(\mathbf{x}_{0},t)$.

######  Lemma 1.

The following asymptotic formula is valid:

(8)
$\displaystyle\mathop{\displaystyle\int}\limits_{S\left(\mathbf{x}_{0},\theta,D\right)}\frac{\text{exp}\left(\text{i}k\sqrt{\left(x-\xi_{1}\right)^{2}+\left(y-\xi_{2}\right)^{2}+\left(z-z_{0}\right)^{2}}\right)}{4\pi\sqrt{\left(x-\xi_{1}\right)^{2}+\left(y-\xi_{2}\right)^{2}+\left(z-z_{0}\right)^{2}}}m\left(\eta,\theta,\boldsymbol{\xi},\mathbf{x}_{0},D\right)d\xi_{1}d\xi_{2}$
$\displaystyle=\begin{cases}\frac{\text{i}}{2k}e^{\text{i}k\left|z-z_{0}\right|}\left[m\left(\eta,\theta,\mathbf{x},\mathbf{x}_{0},D\right)+\mathcal{O}\left(1/k\right)\right]&\text{if
}\left(x,y\right)\in S\left(\mathbf{x}_{0},\theta,D\right),\\\
\mathcal{O}\left(1/k^{2}\right)&\text{if }\left(x,y\right)\notin
S\left(\mathbf{x}_{0},\theta,D\right).\end{cases}$

Proof. Consider polar coordinates for the left-hand side of (8) with the
center at $(x,y)$:

$\displaystyle\xi_{1}-x=r\cos(\varphi),\quad\xi_{2}-y=r\sin(\varphi).$

To simplify the presentation, we take
$m(\xi_{1},\xi_{2})=m\left(\eta,\theta,\boldsymbol{\xi},\mathbf{x}_{0},D\right)$
and denote the left-hand side of (8) by $I_{S}$. Then, we arrive at

$\displaystyle
I_{S}=\int_{0}^{2\pi}\int_{r_{1}\left(\varphi,x,y\right)}^{r_{2}\left(\varphi,x,y\right)}\frac{\exp\left(\text{i}k\sqrt{r^{2}+\left(z-z_{0}\right)^{2}}\right)}{4\pi\sqrt{r^{2}+\left(z-z_{0}\right)^{2}}}m\left(x+r\cos\left(\varphi\right),y+r\sin\left(\varphi\right)\right)rdrd\varphi.$

Notice that

$\displaystyle\left(r_{1}\left(\varphi,x,y\right),\varphi\right)=\begin{cases}\left(0,\varphi\right)&\text{if
}\left(x,y\right)\in\overline{S\left(\mathbf{x}_{0},\theta,D\right)},\\\
\in\partial S\left(\mathbf{x}_{0},\theta,D\right)&\text{if
}\left(x,y\right)\notin\overline{S\left(\mathbf{x}_{0},\theta,D\right)},\end{cases}$
$\displaystyle\left(r_{2}\left(\varphi,x,y\right),\varphi\right)\in\partial
S\left(\mathbf{x}_{0},\theta,D\right).$

Moreover, we have

$\frac{rdr}{\sqrt{r^{2}+\left(z-z_{0}\right)^{2}}}=d\left(\sqrt{r^{2}+\left(z-z_{0}\right)^{2}}\right).$

Therefore, by integration by parts with respect to $r$ twice, we find that

(9) $\displaystyle I_{s}$
$\displaystyle=\frac{1}{\text{i}k}\left[\int_{0}^{2\pi}\frac{1}{4\pi}\exp\left(\text{i}k\sqrt{r_{2}^{2}+\left(z-z_{0}\right)^{2}}\right)m\left(x+r_{2}\cos\left(\varphi\right),y+r_{2}\sin\left(\varphi\right)\right)d\varphi\right]$
$\displaystyle-\frac{1}{\text{i}k}\left[\int_{0}^{2\pi}\frac{1}{4\pi}\exp\left(\text{i}k\sqrt{r_{1}^{2}+\left(z-z_{0}\right)^{2}}\right)m\left(x+r_{1}\cos\left(\varphi\right),y+r_{1}\sin\left(\varphi\right)\right)d\varphi\right]$
$\displaystyle+\mathcal{O}\left(1/k^{2}\right).$

When $(x,y)\notin S\left(\mathbf{x}_{0},\theta,D\right)$, the first two terms
in the right-hand side of (9) vanish due to the definition of $m$ in (4). When
$(x,y)\in S\left(\mathbf{x}_{0},\theta,D\right)$, only the first term vanishes
and thus, we observe that

$\displaystyle-\frac{1}{\text{i}k}\left[\int_{0}^{2\pi}\frac{1}{4\pi}\exp\left(\text{i}k\sqrt{r_{1}^{2}+\left(z-z_{0}\right)^{2}}\right)m\left(x+r_{1}\cos\left(\varphi\right),y+r_{1}\sin\left(\varphi\right)\right)d\varphi\right]$
$\displaystyle=\frac{\text{i}}{2k}e^{\text{i}k\left|z-z_{0}\right|}m\left(x,y\right).$

Thus, (8) follows. $\square$

Lemma 1 implies that for large values of $k$ the incident wave
$v_{0}\left(\mathbf{x},\mathbf{x}_{0},k\right)$ in (7) is indeed basically
contained in the cylinder whose base is the disk
$S\left(\mathbf{x}_{0},\theta,D\right)$. Furthermore, the magnitude of this
field changes with the distance $\left|z-z_{0}\right|$ only as
$\mathcal{O}(1/k)$.

###  3.2 Inverse problem

In our inverse problem, we want to calculate the spatial distribution of the
dielectric constant $\varepsilon_{r}(\mathbf{x})$ in PDE of (2), provided that
(1) holds and our SAR data
$F(\mathbf{x}_{0},t)=u(\mathbf{x}_{0},\mathbf{x}_{0},t)$ are given for
$\mathbf{x}_{0}\in L_{\text{scr}}$, $t>0$. Recall that this is an
underdetermined coefficient inverse problem because the unknown function
depends on three variables, while the data depend only on two variables.
Therefore, our goal of imaging the function $\varepsilon_{r}(\mathbf{x})$ is
extremely limited.

Instead of working on that imaging, we work on finding a certain function
$\widetilde{\varepsilon}_{r}(\mathbf{x})$ on the slant range $P$ (i.e. for
$\mathbf{x}\in P$). It is desirable that this function
$\widetilde{\varepsilon}_{r}(\mathbf{x})$ would characterize well the values
of the actual dielectric constant $\varepsilon_{r}(\mathbf{x})$. Furthermore,
our hope is that some sizes of the computed slant-range image will not be far
from the real sizes of real targets. It is also worth noting that the actual
target is not required to intersect with the slant range $P$. However, it is
still possible to obtain the above parameters of that target because the
radiation patterns of antennas might intersect with this target, thus causing
back reflected waves which arrive at the detectors.

######  Remark 2.

We call $\widetilde{\varepsilon}_{r}(\mathbf{x})$ “the slant range image”.
From the Physics standpoint, the “slant-range image” does not exist in a sense
that it is not a physical object. In Physics, the scattering can be produced
only by the real physical object. The mathematical inversion procedure
described in the present paper takes into account backscattering waves
arriving to any detector along different directions. We use the assumption
that the radiation pattern of the antenna has the same intensity in the slant-
range plane as in the case when the antenna is pointed at the target.

##  4 Imaging via the solution of a coefficient inverse problem for a 1D
hyperbolic PDE

In this part, we propose a way to form the image in the slant range $P$. Prior
to that, we want to mention that the original SAR data
$F(\mathbf{x}_{0},t)=u(\mathbf{x}_{0},\mathbf{x}_{0},t)$ for
$\mathbf{x}_{0}\in L_{\text{scr}}$, $t>0$ need to be preprocessed due to radar
signals. We will discuss our choice of preprocessing procedure in the
numerical section.

Now, let $\mathbf{x}_{0,n}\in L_{\text{scr}}$ be the source location number
$n=\overline{1,N}$. Let $f_{0,n}(t)$ be the preprocessed data of
$F(\mathbf{x}_{0},t)$ at $\mathbf{x}_{0,n}$ for $n=\overline{1,N}$. Consider a
“pseudo” variable $x\in\mathbb{R}$ and let the function $c(x)\in
C^{3}(\mathbb{R})$. We scale to $x\in(0,1)$ the length of the interval of our
interest, i.e. the length of the interval in the slant range of the central
line $CRL(\mathbf{x}_{0,n})$ along which we want to calculate our function
$\widetilde{\varepsilon}_{r}(\mathbf{x})$. Our idea below is to calculate the
slant-range function $\widetilde{\varepsilon}_{r}(\mathbf{x})$ mentioned in
section 3.2 via finding $N$ functions $c(x):=c_{n}(x)$ along each central line
$CRL(\mathbf{x}_{0,n})$ and then forming the 2D image in the slant range plane
out of these $N$ one-dimensional functions $c_{n}(x)$. This idea was initiated
in [21], whereas the related coefficient inverse problem for a 1D hyperbolic
PDE was considered in [37].

###  4.1 Statement of the coefficient inverse problem for a 1D hyperbolic PDE

Given a number $\overline{c}>1$, we assume that

(10) $C^{3}(\mathbb{R})\ni
c\left(x\right)=\begin{cases}\in\left[1,\overline{c}\right]&\text{for
}x\in\left(0,1\right),\\\ 1&\text{for }x\leq 0\text{ and }x\geq 1.\end{cases}$

Let $T>0$. In this scenario, our forward problem is formed as the following
initial value (or Cauchy) problem:

(11) $\begin{cases}c\left(x\right)u_{tt}=u_{xx}&\text{for
}x\in\mathbb{R},t\in\left(0,T\right),\\\ u\left(x,0\right)=0,\quad
u_{t}\left(x,0\right)=\delta\left(x\right)&\text{for
}x\in\mathbb{R}.\end{cases}$

Now we state our targeted inverse problem.

Coefficient Inverse Problem (CIP). _Determine the coefficient_
$c\left(x\right)$ _for_ $x\in\left(0,1\right)$_in the PDE of ( 11), provided
that the following two functions _$f_{0}\left(t\right)$ _and_
$f_{1}\left(t\right)$_are given:_

(12) $u\left(0,t\right)=f_{0}\left(t\right),\quad
u_{x}\left(0,t\right)=f_{1}\left(t\right),\quad t\in\left(0,T\right).$

######  Remark 3.

The data at the point $x=0$ in (12) corresponds to those at the source
position $\mathbf{x}_{0,n}$ in our SAR setting for each $n=\overline{1,N}$.
The Neumann-type data $f_{1}(t)$ in (12) is obtainable. Indeed, it was shown
in [37] (and also commenced in [38]) that the function $u(x,t)$ in the PDE
(11) satisfies the so-called absorbing boundary conditions:

(13) $\displaystyle
u_{x}\left(a_{1},t\right)-u_{t}\left(a_{1},t\right)=0,\quad a_{1}\leq 0,$ (14)
$\displaystyle u_{x}\left(a_{2},t\right)+u_{t}\left(a_{2},t\right)=0,\quad
a_{2}\geq\sqrt{\overline{c}}.$

In (13), take $a_{1}=0$. Then, using (12), we obtain
$f_{1}(t)=f_{0}^{\prime}(t)$. Thus, for each $n=\overline{1,N}$ we find the
Neumann data $f_{1}(t)$ for our CIP by differentiating the preprocessed data
$f_{0,n}(t)$.

######  Remark 4.

After solving the aforementioned CIP, we scale back the “pseudo” interval
$x\in(0,1)$ and obtain $x_{n}^{\prime}\in(0,a)$ for a certain number $a>0$. In
fact, for each $n=\overline{1,N}$ we assign

$\left.\widetilde{\varepsilon}_{r}\left(\mathbf{x}\right)\right|_{\mathbf{x}\in
CRL\left(\mathbf{x}_{0,n}\right)}:=\widetilde{\varepsilon}_{r,n}\left(\mathbf{x}\right):=c\left(x_{n}^{\prime}\right),\quad
x_{n}^{\prime}\in\left(0,a\right).$

We define the function $\widetilde{\varepsilon}_{r}\left(\mathbf{x}\right)$ as
a set of values of the function $\widetilde{\varepsilon}_{r}(\mathbf{x})$ on
1D cross-sections of the slant range $P$ by central lines
$CRL(\mathbf{x}_{0,n})$. Thus, our strategy is to form the 2D image on the
slant range $P$ as the collection of 1D cross-sections:

$\widetilde{\varepsilon}_{r}\left(\mathbf{x}\right)=\left\\{\widetilde{\varepsilon}_{r,n}\left(\mathbf{x}\right),n=\overline{1,N}\right\\}.$

This finalizes our steps of the SAR imaging procedure.

Now we want to discuss the analysis of our CIP. When doing so, we revisit a
transformation investigated in [37] to delineate an essential nonlocal
nonlinear PDE for our mathematical analysis.

Cf. [32], consider the one-to-one change of variables

$y=y(x)=\int_{0}^{x}\sqrt{c(s)}ds.$

Then $y(x)$ is called the travel time, which the wave needs to travel from the
source $\left\\{0\right\\}$ to a certain point $\left\\{x\right\\}$. By (10),
$dy/dx=\sqrt{c(y)}\geq 1$. Moreover, $y(x)=x$ for $x\leq 0$ and $x(y)=y$ for
$y\leq 0$. We introduce the following functions $w,Q,p$:

(15) $\displaystyle
w\left(y,t\right)=u\left(x\left(y\right),t\right)c^{1/4}\left(x\left(y\right)\right),$
(16) $\displaystyle Q\left(y\right)=c^{-1/4}\left(x\left(y\right)\right),\quad
p\left(y\right)=\frac{Q^{\prime\prime}\left(y\right)}{Q\left(y\right)}-2\left[\frac{Q^{\prime}\left(y\right)}{Q\left(y\right)}\right]^{2}.$

Due to (10), functions $Q$ and $p$ are well-defined and we have

$\displaystyle w(y,t)=u(x(y),t)\quad\text{for }y<0\;\text{and
}y>\sqrt{\overline{c}},$ (17) $\displaystyle p\left(y\right)\in
C^{1}\left(\mathbb{R}\right),\quad p\left(y\right)=0\quad\text{for
}y<0\;\text{and }y>\sqrt{\overline{c}}.$

Therefore, we obtain from (11)–(12), (15), (16) that

(18) $\displaystyle w_{tt}=w_{yy}+p\left(y\right)w\quad\text{for
}y\in\mathbb{R},t\in\left(0,\widetilde{T}\right),$ (19) $\displaystyle
w\left(y,0\right)=0,\quad
w_{t}\left(y,0\right)=\delta\left(y\right)\quad\text{for }y\in\mathbb{R},$
(20) $\displaystyle w\left(0,t\right)=f_{0}\left(t\right),\quad
w_{y}\left(0,t\right)=f_{1}\left(t\right)\quad\text{for
}t\in\left(0,\widetilde{T}\right),$

where the number $\widetilde{T}\geq 2\sqrt{\overline{c}}$ depends on $T$.

######  Remark 5.

It was shown in [32, Chapter 2] and [36, Section 2] that the forward problem
(18)–(19) has a unique solution in $C^{3}(\mathbb{R})$. In particular, it
holds that

(21) $\displaystyle
w\left(y,t\right)=\frac{H\left(t-\left|y\right|\right)}{2}+\frac{1}{2}\int_{\left(y-t\right)/2}^{\left(y+t\right)/2}p\left(\xi\right)\left(\int_{\left|\xi\right|}^{t-\left|y-\xi\right|}w\left(\xi,\tau\right)d\tau\right)d\xi,$
(22) $\displaystyle w\left(y,t\right)=0\quad\text{for
}t<\left|y\right|,\quad\lim_{t\rightarrow\left|y\right|^{+}}w\left(y,t\right)=\frac{1}{2},$

where $H(z)$ for $z\in\mathbb{R}$ is the Heaviside function. Here, $w\in
C^{3}(t\geq\left|y\right|)$ is guaranteed by $p(y)\in C^{1}(\mathbb{R})$,
which explains the condition $c(x)\in C^{3}(\mathbb{R})$ imposed in (10).

######  Remark 6.

It is straightforward that to find $c(x)\in C^{3}(\mathbb{R})$ stated in our
CIP, we need to determine the function $p(y)\in C^{1}(\mathbb{R}),$ which is
the solution of the CIP (18)–(20). Since $y\leq\sqrt{\overline{c}}$, we find
$p(y)$ for $y\in(0,b)$ for an arbitrary $b\geq\sqrt{\overline{c}}$. Then the
“threshold” $\widetilde{T}$ for the CIP (18)–(20) is taken by a certain
$\widetilde{T}\geq\max\left\\{2b,T\right\\}\geq 2\sqrt{\overline{c}}$ in order
to guarantee the uniqueness of the CIP (18) and (20); cf. [32].

Let the parameter $\beta\in(0,1/2).$ We will choose it later. Consider an
arbitrary number $\mu\in(0,2\beta b)$. We introduce the following rectangle
$R$ and the triangle $R_{\beta,\mu,b}\subset R$, as follows:

(23) $\displaystyle
R:=\left\\{\left(y,t\right)\in\mathbb{R}^{2}:y\in\left(0,b\right),t\in\left(0,\widetilde{T}\right)\right\\},$
(24) $\displaystyle
R_{\beta,\mu,b}:=\left\\{\left(x,t\right)\in\mathbb{R}^{2}:x+\beta t<2\beta
b-\mu\;\text{for }x,t>0\right\\}.$

Henceforth, we consider a new function $v(y,t)$,

(25) $v\left(y,t\right)=w\left(y,t+y\right),\quad x,t>0.$

Then $v\in C^{3}(y\geq 0,t\geq 0)$ and (21)–(22) imply

(26) $\lim_{t\rightarrow 0^{+}}v\left(y,t\right)=\lim_{t\rightarrow
0^{+}}w\left(y,t+y\right)=\lim_{t\rightarrow
y^{+}}w\left(y,t\right)=\frac{1}{2}.$

It follows from (18) and (26) that

(27) $\displaystyle v_{yy}-2v_{yt}+p\left(y\right)v=0\quad\text{for
}\left(y,t\right)\in R,$ (28) $\displaystyle
v\left(y,0\right)=\frac{1}{2}\quad\text{for }y\in\left(0,b\right).$

Take $t=0$ in (27) and use (28), (26). We deduce that

(29) $p(y)=4v_{yt}(y,0),\quad y\in(0,b).$

Next, we differentiate both sides of (27) with respect to $t$ and consider
$V(y,t)=v_{t}(y,t)$ satisfying $V\in C^{2}(\overline{R})$. Using (12), (14)
and (25)–(29), we obtain the following nonlocal nonlinear boundary value
problem for a hyperbolic equation:

(30) $\displaystyle V_{yy}-2V_{yt}+4V_{y}\left(y,0\right)V=0\quad\text{for
}\left(y,t\right)\in R,$ (31) $\displaystyle
V\left(0,t\right)=q_{0}\left(t\right):=f_{0}^{\prime}\left(t\right)\quad\text{for
}t\in\left(0,\widetilde{T}\right),$ (32) $\displaystyle
V_{y}\left(0,t\right)=q_{1}\left(t\right):=f_{0}^{\prime\prime}\left(t\right)+f_{1}^{\prime}\left(t\right)\quad\text{for
}t\in\left(0,\widetilde{T}\right),$ (33) $\displaystyle
V_{y}\left(b,t\right)=0\quad\text{for }t\in\left(0,\widetilde{T}\right).$

If $V(y,t)$ for $(y,t)\in R$ is found via solving the overdetermined boundary
value problem (30)–(33), then we can rely on (29) to have

(34) $p(y)=4V_{y}(y,0),\quad y\in(0,b).$

Then recalling Remark 6 we are able to determine the coefficient $c(x)$ of the
original CIP; see (16). Thus, we need to solve the nonlinear overdetermined
boundary value problem (30)–(33).

###  4.2 Convexification method for the nonlocal nonlinear problem

Consider the rectangle $R$ in (23). Whenever we deal below with function
spaces $H^{k}(R)$, we work only with real-valued functions. Consider now the
subspaces $H_{0}^{2}(R)\subset H^{2}(R)$ and $H_{0}^{4}(R)\subset H^{4}(R)$,

$\displaystyle H_{0}^{2}\left(R\right):=\left\\{u\in
H^{2}\left(R\right):u\left(0,t\right)=u_{y}\left(0,t\right),u_{y}\left(b,t\right)=0\right\\},$
$\displaystyle H_{0}^{4}\left(R\right):=H^{4}\left(R\right)\cap
H_{0}^{2}\left(R\right),$

see Remark 7 about $H^{4}(R).$ Recall that the parameters $\beta\in(0,1/2)$,
$b\geq\sqrt{\overline{c}}$ are introduced in (24). As introduced in [36, 37],
we consider $\psi_{\lambda}(y,t)$,

(35) $\psi_{\lambda}(y,t)=e^{-2\lambda(y+\beta t)},\quad\lambda\geq 1.$

This is the Carleman Weight Function (CWF) for the linear operator
$\partial_{y}^{2}-2\partial_{y}\partial_{t}$. This operator is the principal
part of the operator of equation (30), which we want to solve.

######  Theorem 2 (Carleman estimate [36]).

There exist constants $C=C(\beta,R)>0$ and
$\lambda_{0}=\lambda_{0}(\alpha,R)\geq 1$ depending only on listed parameters
such that for all functions $u\in H_{0}^{2}(R)$ and for all
$\lambda\geq\lambda_{0}$ the following Carleman estimate is valid:

$\displaystyle\int\limits_{R}\left(u_{yy}-2u_{yt}\right)^{2}\psi_{\lambda}dydt\geq
C\int\limits_{R}\left(\lambda\left(u_{y}^{2}+u_{t}^{2}\right)+\lambda^{3}u^{2}\right)\psi_{\lambda}dydt$
$\displaystyle+C\int\limits_{0}^{b}\left(\lambda
u_{y}^{2}+\lambda^{3}u^{2}\right)\left(y,0\right)e^{-2\lambda y}dy-
Ce^{-2\lambda\widetilde{T}}\int\limits_{0}^{b}\left(\lambda
u_{y}^{2}+\lambda^{3}u^{2}\right)\left(y,\widetilde{T}\right)dy.$

Now we introduce a Tikhonov-like cost functional $J_{\lambda,\gamma}$ weighted
by the CWF (35), where $\gamma\in(0,1)$ is the regularization parameter. Let
$M>0$ be an arbitrary number. Consider the convex set $B(M,q_{0},q_{1})\subset
H^{4}\left(R\right)$ of the diameter $2M$,

(36) $\displaystyle B\left(M,q_{0},q_{1}\right)$ $\displaystyle=\left\\{u\in
H^{4}\left(R\right):u\left(0,t\right)=q_{0}\left(t\right),u_{y}\left(0,t\right)=q_{1}\left(t\right),u_{y}\left(b,t\right)=0,\left\|u\right\|_{H^{4}\left(R\right)}<M\right\\}.$

In (36), $q_{0}$ and $q_{1}$ are given data in (31) and (32), respectively.

######  Remark 7.

We need to use the space $H^{4}(R)$ since $H^{4}(R)\subset
C^{2}(\overline{R})$ is known by the embedding theorem. Theorem 3 is a direct
analog of Theorem 2 of [37], which requires the regularity
$C^{2}(\overline{R})$; see the inequality (73) in [37]. However, in
computations we use $H^{2}(R)$ instead of $H^{4}(R)$, see section 5. Below,
$B\left(M,q_{0},q_{1}\right)$ denotes only interior points of the set
$\overline{B\left(M,q_{0},q_{1}\right)}$, which is closed in the space
$H^{4}(R)$.

Following (30), denote

(37) $\mathcal{S}(V)=V_{yy}-2V_{yt}+4V_{y}(y,0)V,\quad(y,t)\in R.$

Henceforth, our weighted Tikhonov-like cost functional
$J_{\lambda,\gamma}:H^{5}(R)\rightarrow\mathbb{R}_{+}$ is defined as follows:

(38)
$J_{\lambda,\gamma}\left(V\right)=\mathop{\displaystyle\int}\limits_{R}\left[\mathcal{S}\left(V\right)\right]^{2}\psi_{\lambda}dydt+\gamma\left\|V\right\|_{H^{4}\left(R\right)}^{2},\quad\gamma\in(0,1).$

Minimization Problem. _Minimize the Tikhonov-like cost functional_
$J_{\lambda,\gamma}\left(V\right)$_in ( 38) on the set
_$\overline{B\left(M,q_{0},q_{1}\right)}.$

######  Theorem 3 (Global strict convexity [37]).

For all $\lambda,\gamma>0$ and for all elements
$V\in\overline{B\left(M,q_{0},q_{1}\right)}$, the cost functional (38) has the
Fréchet derivative $J_{\lambda,\gamma}^{\prime}\left(V\right)\in
H_{0}^{4}\left(R\right).$ Let $\lambda_{0}=\lambda_{0}(\beta,R)\geq 1$ be the
number obtained in Theorem 2. Let $R_{\alpha,\mu,b}$ be the triangle defined
in (24) and let $\widetilde{T}\geq 2b.$ Then there exists a sufficiently large
number $\lambda_{1}=\lambda_{1}\left(R,M,\beta\right)\geq$ $\lambda_{0}$
depending only on $R,M$ and $\beta$ such that for all $\lambda\geq\lambda_{1}$
and for all $\gamma\in[2e^{-\lambda\beta\widetilde{T}},1),$ the functional
$J_{\lambda,\gamma}(V)$ is strictly convex on the set
$\overline{B\left(M,q_{0},q_{1}\right)},$ i.e. for all
$V_{1},V_{2}\in\overline{B\left(M,q_{0},q_{1}\right)}$ the following estimate
holds

(39) $\displaystyle
J_{\lambda,\gamma}\left(V_{2}\right)-J_{\lambda,\gamma}\left(V_{1}\right)-J_{\lambda,\gamma}^{\prime}\left(V_{1}\right)\left(V_{2}-V_{1}\right)\geq
Ke^{-2\lambda\left(2\beta
b-\mu\right)}\left\|V_{2}-V_{1}\right\|_{H^{1}\left(R_{\beta,\mu,b}\right)}^{2}$
$\displaystyle+Ke^{-2\lambda\left(2\beta
b-\mu\right)}\left\|V_{2}\left(y,0\right)-V_{1}\left(y,0\right)\right\|_{H^{1}\left(0,2\beta
b-\mu\right)}^{2}+\frac{\gamma}{2}\left\|V_{2}-V_{1}\right\|_{H^{4}\left(R\right)}^{2},$

where the constant $K=K\left(\beta,\mu,b,M\right)>0$ depends only on listed
parameters.

To this end, $K>0$ denotes different constants depending on some parameters,
to be specified. Theorem 3 claims the global strict convexity since smallness
conditions are not imposed on the diameter $2M>0$ of the convex set
$B(M,q_{0},q_{1})$. Below we formulate the theorem about the existence and
uniqueness of a minimizer of the functional $J_{\lambda,\gamma}(V)$ on the set
$\overline{B(M,q_{0},q_{1})}$.

######  Theorem 4 (Existence and uniqueness of a minimizer [37]).

Let parameters $\lambda_{1}$ and $\gamma$ be the same as in Theorem 3. Then
for any $\lambda\geq\lambda_{1}$ there exists a unique minimizer
$V_{\min,\lambda,\gamma}$ of the functional $J_{\lambda,\gamma}\left(V\right)$
on the set $\overline{B\left(M,q_{0},q_{1}\right)}.$ Furthermore, the
following inequality holds

(40)
$J_{\lambda,\gamma}^{\prime}\left(V_{\min,\lambda,\gamma}\right)\left(V-V_{\min,\lambda,\gamma}\right)\geq
0\quad\text{for any }V\in\overline{B\left(M,q_{0},q_{1}\right)}.$

Using the standard concept of the regularization theory (cf. [6, 40]), we
assume the existence of ideal, noiseless data
$f_{0}^{\ast}\left(t\right),f_{1}^{\ast}\left(t\right)$ in (32), which
correspond to the exact coefficient $p^{\ast}\left(y\right)\in
C^{1}\left(\mathbb{R}\right)$ and the exact function
$w^{\ast}\left(y,t\right)$ in (18)–(20). In addition, we assume that the
function $p^{\ast}\left(y\right)$ satisfies conditions (17). Let the function
$V^{\ast}\left(y,t\right)$ be obtained from the function
$w^{\ast}\left(y,t\right)$ by the same way we have obtained above $V$ from
$w$; cf. section 4.1. Assume that

(41) $\displaystyle V^{\ast}\in B(M,q_{0}^{\ast},q_{1}^{\ast}),$

where functions $q_{0}^{\ast},q_{1}^{\ast}$ are the noiseless data of
$q_{0},q_{1}$ in (31), (32), respectively. By (30)–(33) and (37), we have

$\displaystyle\mathcal{S}(V^{\ast})=0,\quad(y,t)\in R,\;\text{and
}p^{\ast}(y)=4v_{y}^{\ast}(y,0),\quad y\in(0,b).$

Let $\sigma>0$ be a sufficiently small number which we will choose later. This
number characterizes the level of noise in the data $q_{0},q_{1}.$ Assume that
$\sigma\in\left(0,\min\left\\{1,M\right\\}\right)$. Suppose that there exist
functions $F\in B\left(M,q_{0},q_{1}\right)$ and $F^{\ast}\in
B\left(M,q_{0}^{\ast},q_{1}^{\ast}\right)$ such that

(42) $\left\|F-F^{\ast}\right\|_{H^{4}(R)}<\sigma.$

The presence of functions $F$ and $F^{\ast}$ is to obtain the zero Dirichlet
and Neumann boundary condition at $\left\\{x=0\right\\}$ in (31) and (32). In
fact, for every $V\in B\left(M,q_{0},q_{1}\right)$ we consider

$W^{\ast}=V^{\ast}-F^{\ast},\quad W=V-F.$

Then by the triangle inequality, it is clear that

$W,W^{\ast}\in B_{0}\left(2M\right):=\left\\{W\in
H_{0}^{4}\left(R\right):\left\|W\right\|_{H^{4}\left(R\right)}<2M\right\\}.$

Moreover, for any $W\in B_{0}\left(2M\right)$ it holds that $W+F\in
B\left(3M,q_{0},q_{1}\right)$. Therefore, we can consider a modification
$\tilde{J}_{\lambda,\gamma}:B_{0}(2M)\rightarrow\mathbb{R}_{+}$ of the
functional $J_{\lambda,\gamma}$,

$\tilde{J}_{\lambda,\gamma}(W)=J_{\lambda,\gamma}(W+F)\quad\text{for any }W\in
B_{0}(2M).$

As inherited from properties of $J_{\lambda,\gamma}$, the functional
$\tilde{J}_{\lambda,\gamma}$ is globally strict convex on the ball
$\overline{B_{0}(2M)}\subset H_{0}^{4}(R)$ and admits a unique minimizer
$W_{\min,\lambda,\gamma}$ on that ball. For brevity, we do not formulate the
theorem for $\tilde{J}_{\lambda,\gamma}$ as it is analogous to the statement
in Theorem 3. The presence of the functional $\tilde{J}_{\lambda,\gamma}$ is
necessary to estimate the accuracy of the minimizer $V_{\min,\lambda,\gamma}$
in the next theorem; see [37, Theorem 5] for details of the proof.

######  Theorem 5 (Accuracy estimate of the minimizer [37]).

Suppose that conditions (41) and (42) hold. Assume that $\widetilde{T}\geq
4b.$ Denote

$\omega=\frac{\alpha\left(\widetilde{T}-4b\right)+\mu}{2\left(2\beta
b-\mu\right)},\quad\rho=\frac{1}{2}\min\left\\{\omega,1\right\\}.$

Let $\lambda_{1}=\lambda_{1}\left(R,M,\beta\right)$ be the number of Theorem
3. Denote $\lambda_{2}=\lambda_{1}\left(R,3M,\beta\right)\geq\lambda_{1}.$ Let
the number $\sigma_{0}\in\left(0,1\right)$ be sufficiently small such that
$\ln\sigma_{0}^{-1/\left(2\left(2\beta b-\mu\right)\right)}\geq\lambda_{2}.$
Let the number $\sigma\in\left(0,\sigma_{0}\right).$ Let the numbers
$\lambda=\lambda\left(\sigma\right)$ and $\gamma=\gamma\left(\sigma\right)$ be
defined as:

(43)
$\displaystyle\lambda=\lambda\left(\sigma\right)=\ln\sigma^{-1/\left(2\left(2\beta
b-\mu\right)\right)}>\lambda_{2},$ (44)
$\displaystyle\gamma=\gamma\left(\sigma\right)=2e^{-\lambda\beta\widetilde{T}}=2\sigma^{\left(\beta\widetilde{T}\right)/\left(2\left(2\beta
b-\mu\right)\right)}.$

Let $V_{\min,\lambda,\gamma}\left(y,t\right)$ be the unique minimizer of the
cost functional $J_{\lambda,\gamma}\left(V\right)$ on the set
$B\left(M,q_{0},q_{1}\right)$, cf. Theorem 4. Let
$p_{\min,\lambda,\gamma}\left(y\right)$ be the coefficient obtained from the
minimizer $V_{\min,\lambda,\gamma}\left(y,t\right)$ via (34), i.e.
$p_{\min,\lambda,\gamma}\left(y\right)=4\partial_{y}V_{\min,\lambda,\gamma}\left(y,0\right)$.
Then the following accuracy estimates are valid

(45)
$\displaystyle\left\|V_{\min,\lambda,\gamma}-V^{\ast}\right\|_{H^{1}\left(R_{\beta,\mu,b}\right)}\leq
K\sigma^{\rho},$ (46)
$\displaystyle\left\|p_{\min,\lambda,\gamma}-p^{\ast}\right\|_{L^{2}\left(0,2\beta
b-\mu\right)}\leq K\sigma^{\rho}.$

######  Theorem 6 (Lipschitz continuity of $J_{\lambda,\gamma}^{\prime}$).

As in Theorem 3, let $J_{\lambda,\gamma}^{\prime}\left(v\right)\in
H_{0}^{4}\left(R\right)$ be the Fréchet derivative of the functional
$J_{\lambda,\gamma}\left(V\right)$. Then for any $\lambda,\gamma>0$ it
satisfies the Lipschitz continuity condition on the set
$\overline{B\left(M,q_{0},q_{1}\right)}.$ More precisely, there exists a
number $\widetilde{C}=\widetilde{C}\left(\lambda,\gamma,M\right)>0$ depending
only on listed parameters such that

$\left\|J_{\lambda,\gamma}^{\prime}\left(V_{1}\right)-J_{\lambda,\gamma}^{\prime}\left(V_{2}\right)\right\|_{H^{4}\left(R\right)}\leq\widetilde{C}\left\|V_{1}-V_{2}\right\|_{H^{4}\left(R\right)}\text{
for all }V_{1},V_{2}\in\overline{B\left(M,q_{0},q_{1}\right)}.$

The proof of this theorem is omitted here because it is completely similar
with the proof of Theorem 3.1 of [5].

######  Remark 8.

In our first paper on SAR imaging [21], we have seen a large amount of
numerical and experimental results using the gradient descent method, although
the convergence theory is based on the gradient projection method in [37]. The
only reason is that the gradient descent method is rather easy-to-implement
numerically. Therefore, this is the moment to explore the convergence analysis
of the gradient descent method. It is also accentuated that starting from the
work [5], only theorems about the global convergence of the gradient
projection method in the context of the convexification method were studied;
cf. e.g. [14, 15, 16, 18, 20, 36, 37]. Hence, our analysis below is really
important to fulfill the gap between the theory and numerics of the
convexification method.

###  4.3 Global convergence of the gradient descent method

Let the number $\varkappa\in(0,1)$ and the number $\vartheta\in(0,1/3)$. Let
$V^{(0)}\in B\left(\vartheta M,q_{0},q_{1}\right)$ be an arbitrary point in
that set. We set the gradient descent method as follows:

(47) $V^{\left(n\right)}=V^{\left(n-1\right)}-\varkappa
J_{\lambda,\gamma}^{\prime}\left(V^{\left(n-1\right)}\right),\quad
n=1,2,\ldots$

Note that $V^{(n)}$ denotes the $n$th iteration that expects to be close to
the minimizer $V_{\min,\lambda,\gamma}$ as $n$ tends to infinity. It is also
worth noting that since $J_{\lambda,\gamma}^{\prime}\left(V^{(n-1)}\right)\in
H_{0}^{4}(R)$ by Theorem 3, then all functions $V^{(n)}$ of the sequence (47)
satisfy the same boundary conditions as ones in (36).

As assumed in (41), we know that the exact solution $V^{\ast}$ is an interior
point of the set $B(M,q_{0}^{\ast},q_{1}^{\ast})$. This is an interior point
in terms of the norm of the space $H^{4}(R)$. On the other hand, by (45) the
distance between the exact solution $V^{\ast}$ and the minimizer
$V_{\min,\lambda,\gamma}$ is small in the norm of the space
$H^{1}\left(R_{\beta,\mu,b}\right)$, as long as the noise level $\sigma$ in
the data is sufficiently small. Therefore, it is reasonable to assume that the
minimizer $V_{\min,\lambda,\gamma}$ is an interior point of the set
$B(M,q_{0},q_{1})$ rather than the one located on the boundary of this set,
see Theorem 4. More precisely, we need to assume that

(48) $V_{\min,\lambda,\gamma}\in B\left(\vartheta M,q_{0},q_{1}\right).$

Now we are in position to estimate the distance between the approximation
$V^{(n)}$ and the minimizer $V_{\min,\lambda,\gamma}$.

######  Theorem 7 (Global convergence of the gradient method).

Fix a number $\vartheta\in(0,1/3)$. Assume that (48) holds with this
$\vartheta.$ Suppose that the parameters $\lambda$ and $\gamma$ are the same
as in Theorem 3. Then there exists a sufficiently small number
$\varkappa_{0}\in(0,1)$ such that the sequence
$\left\\{V^{(n)}\right\\}_{n=0}^{\infty}\subset B\left(M,q_{0},q_{1}\right)$
for each $\varkappa\in(0,\varkappa_{0})$. Furthermore, for each
$\varkappa\in(0,\varkappa_{0})$ there exists a number
$\theta=\theta(\varkappa)\in(0,1)$ such that

(49)
$\left\|V_{\min,\lambda,\gamma}-V^{\left(n\right)}\right\|_{H^{4}\left(R\right)}\leq\theta^{n}\left\|V_{\min,\lambda,\gamma}-V^{\left(0\right)}\right\|_{H^{4}\left(R\right)},\quad
n=1,2,\ldots$

Proof. Consider the operator $Z:B(M,q_{0},q_{1})\rightarrow H^{4}(R)$ defined
as

(50) $Z(V)=V-\varkappa J_{\lambda,\gamma}^{\prime}(V).$

It follows from the proof of [5, Theorem 2.1] that Theorems 3 and 6 imply that
there indeed exist numbers $\varkappa_{0}\in(0,1)$ and
$\theta(\varkappa)\in(0,1)$ for $\varkappa\in(0,\varkappa_{0})$ such that the
operator $Z$ defined in (50) has the contractual mapping property, which reads
as

(51)
$\left\|Z\left(V_{1}\right)-Z\left(V_{2}\right)\right\|_{H^{4}\left(R\right)}\leq\theta\left\|V_{1}-V_{2}\right\|_{H^{4}\left(R\right)}\text{
for all }V_{1},V_{2}\in B\left(M,q_{0},q_{1}\right).$

Since $B\left(\vartheta M,q_{0},q_{1}\right)\subset
B\left(M,q_{0},q_{1}\right)$, then estimate (51) holds for all $V_{1},V_{2}\in
B\left(\vartheta M,q_{0},q_{1}\right)$. However, it was not proven in [5] that
the operator $Z$ maps the set $B\left(M,q_{0},q_{1}\right)$ into itself. Note
that since $J_{\lambda,\gamma}^{\prime}\in H_{0}^{4}(R)$, then functions
$Z(V)$ and $V$ satisfy the same boundary conditions for any $V\in
B\left(M,q_{0},q_{1}\right)$.

Since $V_{\min,\lambda,\gamma}$ is an interior point of the set
$B\left(\vartheta M,q_{0},q_{1}\right),$ then
$J_{\lambda,\gamma}^{\prime}\left(V_{\min,\lambda,\gamma}\right)=0$.
Therefore,

$V_{\min,\lambda,\gamma}=V_{\min,\lambda,\gamma}-\varkappa
J_{\lambda,\gamma}^{\prime}\left(V_{\min,\lambda,\gamma}\right)=Z\left(V_{\min,\lambda,\gamma}\right)\in
B\left(\vartheta M,q_{0},q_{1}\right).$

It now remains to check that $V^{(n)}\in B\left(M,q_{0},q_{1}\right)$ whenever
$V^{(0)}\in B\left(\vartheta M,q_{0},q_{1}\right)$ for any
$\vartheta\in(0,1/3)$. In fact, it can be done by mathematical induction. For
$n=1$, we particularly have

(52) $\displaystyle V^{(1)}=V^{(0)}-\varkappa
J_{\lambda,\gamma}^{\prime}(V^{(0)})=Z(V^{(0)}).$

Therefore, using (51), we estimate that

(53)
$\displaystyle\left\|V^{\left(1\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}\leq\theta\left\|V^{\left(0\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}.$

Recall that $V^{(0)}\in B\left(\vartheta M,q_{0},q_{1}\right)$. By the
triangle inequality, we obtain

$\displaystyle\left\|V^{\left(1\right)}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq\left\|V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}+\theta\left\|V^{\left(0\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq\left(1+\theta\right)\left\|V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}+\theta\left\|V^{\left(0\right)}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq 2\vartheta M+\theta\vartheta M<3\vartheta M<M.\text{ }$

Henceforth, $V^{(1)}\in B(M,q_{0},q_{1})$. In the same vein as (52), we have
for $n=2$ that $V^{(2)}=Z(V^{(1)})$. Therefore, using (51) and (53), we
estimate that

(54)
$\displaystyle\left\|V^{\left(2\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}\leq\theta\left\|V^{\left(1\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}\leq\theta^{2}\left\|V^{\left(0\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}.$

Applying again the triangle inequality and using (54), we obtain

$\displaystyle\left\|V^{\left(2\right)}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq\left\|V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}+\theta^{2}\left\|V^{\left(0\right)}-V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq\left(1+\theta^{2}\right)\left\|V_{\min,\lambda,\gamma}\right\|_{H^{4}\left(R\right)}+\theta^{2}\left\|V^{\left(0\right)}\right\|_{H^{4}\left(R\right)}$
$\displaystyle\leq 2\vartheta M+\theta^{2}\vartheta M<3\vartheta M<M.\text{ }$

Henceforth, $V^{(2)}\in B(M,q_{0},q_{1})$.

By continuing this process up to the $n$-th iteration, we obtain the same that
$V^{(n)}\in B(M,q_{0},q_{1})$ and (49) holds true. Hence, we complete the
proof of the theorem. $\square$

In the following, we prove that the sequence (47) of the gradient descent
method actually converges to the exact solution, as long as the noise level
$\sigma$ in the data (cf. (42)) tends to zero. Then the convergence of the
sequence $\left\\{p^{(n)}\right\\}_{n=0}^{\infty}$, where
$p^{(n)}(y)=4\partial_{y}V_{y}^{(n)}(y,0)$ (see (34)), follows immediately.

######  Theorem 8.

Suppose that conditions of Theorems 5 and 7 hold true. Then the following
convergence estimate is valid:

$\displaystyle\left\|V^{\ast}-V^{\left(n\right)}\right\|_{H^{1}\left(R_{\alpha,\mu,b}\right)}$
$\displaystyle+\left\|p^{\ast}-p^{\left(n\right)}\right\|_{L^{2}\left(0,2\alpha
b-\mu\right)}$ $\displaystyle\leq
2K\sigma^{\rho}+4C\theta^{n}\left\|V_{\min,\lambda,\gamma}-V^{\left(0\right)}\right\|_{H^{4}\left(R\right)},\quad
n=1,2,\ldots,$

where $C=C(R)>0$ depends only on the rectangle $R$.

Proof. Combine (45) and (49). Since $H^{4}(R)\subset
H^{1}\left(R_{\beta,\mu,b}\right)$, then we use the triangle inequality,

$\displaystyle\left\|V^{\ast}-V^{\left(n\right)}\right\|_{H^{1}\left(R_{\beta,\mu,b}\right)}$
$\displaystyle\leq\left\|V^{\ast}-V_{\min,\lambda,\gamma}\right\|_{H^{1}\left(R_{\beta,\mu,b}\right)}+\left\|V_{\min,\lambda,\gamma}-V^{\left(n\right)}\right\|_{H^{1}\left(R_{\beta,\mu,b}\right)}$
$\displaystyle\leq
K\sigma^{\rho}+\theta^{n}\left\|V_{\min,\lambda,\gamma}-V^{\left(0\right)}\right\|_{H^{4}\left(R\right)}.$

We estimate the distance between $p^{\ast}$ and $p^{(n)}$ similarly. Indeed,
we have

$\displaystyle\left\|p^{\ast}-p^{\left(n\right)}\right\|_{L^{2}\left(0,2\beta
b-\mu\right)}$
$\displaystyle\leq\left\|p^{\ast}-p_{\min,\lambda,\gamma}\right\|_{L^{2}\left(0,2\beta
b-\mu\right)}+\left\|p_{\min,\lambda,\gamma}-p^{\left(n\right)}\right\|_{L^{2}\left(0,2\beta
b-\mu\right)}$ $\displaystyle\leq
K\sigma^{\rho}+4\left\|\partial_{y}V_{\min,\lambda,\gamma}\left(\cdot,0\right)-\partial_{y}V^{\left(n\right)}\left(\cdot,0\right)\right\|_{L^{2}\left(0,2\alpha
b-\mu\right)}$ $\displaystyle\leq
K\sigma^{\rho}+4C\theta^{n}\left\|V_{\min,\lambda,\gamma}-V^{\left(0\right)}\right\|_{H^{4}\left(\Omega\right)},$

where we have used (46), (49) and the trace inequality

$\left\|\mathcal{U}_{y}\left(\cdot,0\right)\right\|_{L^{2}\left(0,2b-\mu\right)}\leq
C\left\|\mathcal{U}\right\|_{H^{4}\left(\Omega\right)}\text{ for any
}\mathcal{U}\in H^{4}(R).$

This completes the proof of the theorem. $\square$

######  Remark 9.

Theorem 7 is our central theorem since this is the first time we prove the
global convergence of the gradient descent method by restricting the minimizer
in a smaller set of $B\left(M,q_{0},q_{1}\right)$. It is worth mentioning that
the convergence in Theorem 7 as well as its consequence in Theorem 8 are still
global since the starting point $V^{(0)}$ is arbitrary in $B\left(\vartheta
M,q_{0},q_{1}\right)$ with $\vartheta$ being fixed in $(0,1/3)$, and smallness
conditions are not imposed on the number $M$.

##  5 Numerical verification using experimental data

In our first paper on SAR imaging [21], we verified the proposed method using
the simulated data. This means the SAR data were computed by solving the
Lippmann–Schwinger equation (7) obtained after numerous assumptions about the
forward problem. This equation can be solved using a trigonometric Galerkin-
based numerical approach postulated in [24]. We also tested in [21] our method
using experimental data for a two-sided scanning of a barrack building in
which experimental objects with high dielectric constants were located.

In this work, we want to corroborate the method using experimental data
collected at the University of North Carolina at Charlotte. We have placed
horizontally the watered bottle inside of a wooden framed box fully filled
with the dry sand. The front and back sides of that sandbox were covered by
styrofoam layers. The thickness of each layer is about 7 centimeters (cm),
while the sandbox itself has the dimensions of approximately $50\times
50\times 70$ $\text{cm}^{3}$ (H$\times$W$\times$L). That styrofoam has almost
the same dielectric constant as the air and thus, it does not affect our
backscattering signals. Following the purpose of the search of anti-personnel
landmines, the bottle is buried close to the front surface with 10 cm depth.
Getting a landmine with the 10 cm depth is usually considered in many models
of de-mining operations; cf. e.g. [11].

Figure 2: A photograph of our experimental configuration with a fake landmine
buried in a sandbox. That sandbox was moved along a straight line which has
formed the $\pi/4$ angle with the wall which can be seen in the photograph.
The total length of this movement was 140 cm with the 5 cm step size. We have
used the concurrent movement along a straight line of both the antenna and
detector, see Figure 1. This models well the SAR data collection of Figure 1.

The fake landmine under consideration was already tested in [14, 15] for a
different mathematical model in the frequency domain. Also different from
[14], this time we conduct our experiment in a room that minimizes the
presence of other devices and items. This is reasonable for both de-mining
operations in the deserts and detection services for potential threats using
special drones. However, it can be seen in Figure 2 that we still have many
unwanted obstacles and furniture that can affect the quality of the
backscattering data which we measure. Mounted in front of the sandbox is a
standard horn antenna with 20 cm in length, and our detector is a point probe
placed behind that transmitter. Besides, the distance between the antenna and
the receiver is about 17 cm. In order to obtain the elevation angle, the
sandbox is placed and moved diagonally to the direction that the antenna runs
along with; see Figures 2 and 3. Note that due to technical difficulties
during the experiment, we actually move the sandbox instead of running the
antenna along a straight line. It still obeys the model of the stripmap SAR
imaging introduced in Figure 1 with a particular elevation angle of $\pi/4$.

In our experiment, the detector is moved in a square of $20\times 20$
($\text{cm}^{2}$) with the step size of 2 cm for each source position. This is
because we want the transmitters and the detectors to stay (relatively) in the
same position in the aircraft. This circumstance is different from what we had
in [14, 15], where the measurement square is of $1\times 1$ in $\text{m}^{2}$.
Nonetheless, such a smaller working area of the detector would hinder us from
collecting good-quality signals, since the antenna, unfortunately, is much
larger than the point probe. Among many distractions, this is the significant
challenge when we work with the experimental data in this SAR setting.

(a)  (b)

Figure 3: (a) A schematic diagram for our experimental setup to collect the
SAR data, which are backscattered by a buried target. (b) Basic geometry of
the experimental data collection. In the measurement plane, the red star
denotes the optimal detector determined in the data preprocessing procedure.

It is worth noting that our collected raw data are frequency-dependent, while
SAR imaging works with the time-dependent ones. This circumstance is different
from what we have done in [14, 15, 25], where we have tackled with the
frequency domain data using the Helmholtz equation. Therefore, we apply the
inverse Fourier transform to get the time-dependent data, which can be
understood via the following operation. Recall from section 4 that the
original SAR data is denoted by $F(\mathbf{x}_{0},t)$. We then denote
correspondingly the raw frequency-dependent data by
$\tilde{F}(\mathbf{x}_{0},k)$, where $k>0$ is the wavenumber. The inverse
Fourier transform reads as

(55)
$F(\mathbf{x}_{0},t)=\mathop{\displaystyle\int}\limits_{-\infty}^{\infty}\tilde{F}(\mathbf{x}_{0},k)e^{-\text{i}kt}dk.$

In practice, we actually integrate in (55) over the interval $(k_{1},k_{2})$
defined in section 3.1.

To this end, “m” stands for meter; “s” and “ns” mean second and nanosecond,
respectively. We remark that the relation between the wavenumber $k$ and the
frequency (denoted by Fre) is $k=2\pi\text{Fre}/c_{0}$, where $c_{0}=3\times
10^{8}$ (m/s) is the speed of light in vacuum. In this regard, we approximate
the inverse Fourier transform using the frequencies ranging from 5.6 GHz to
9.0 GHz.

As in [21], we compute an approximation to the minimizer of
$J_{\lambda,\gamma}$ in (38) using the standard finite difference setting. In
this sense, the finite difference operators are applied to approximate the
corresponding differential operators involved in $\mathcal{S}$; see (38). For
brevity, we do not introduce those, but refer to, e.g., [14] and references
cited therein, where we have introduced analogs of discrete operators. Then
the discrete solution is sought by the minimization of the corresponding
discrete cost functional. As in [21], we only consider the simpler penalty
term in $H^{2}(R)$ instead of $H^{4}(R)$ presented in the theoretical section.
This replacement helps to reduce the complexity of computations.

It is worth mentioning that in the same vein as what we have done in [21], we
scale the simulated data by a calibration factor $CF$. This is because the
problem under consideration is a quite challenging one that requires the
settings in the forward and inverse problems to be chosen differently. We hope
to address this in the future.

Consider now the speed of light $c_{0}=0.3$ (m/ns). In the case of dimensions,
the wave equation has the following form:

(56) $\frac{c(y)}{c_{0}^{2}}u_{ss}=u_{yy},$

where we consider $s$ in ns and $y$ in m. Then let
$x=y/\left(0.3\text{m}\right)$ and $t=s/\text{ns}$. We find that

$u_{xx}=\left(0.3\text{m}\right)^{2}u_{yy},\quad
u_{tt}=\left(1\text{ns}\right)^{2}u_{ss}.$

Henceforth, we arrive at a new form of (56) in the dimensionless regime:

(57) $\tilde{c}(x)u_{tt}=u_{xx},\quad\tilde{c}(x)=c(0.3x).$

This equation resembles (11). In this sense, $x=1$ means $y=0.3$ m in the
reality, and $t=1$ means $s=1$ ns. Thus, to construct the image from
experimental data, we solve via the convexification the CIP for equation (57)
with the initial conditions as in the second line of (11), the data for the
CIP as in (12) and the unknown coefficient $\tilde{c}(x).$

In our experiment, the measured data are actually far-field ones. Recall that
outside of the computational domain, which is our sandbox’s length here, the
dielectric constant equals to unity. Therefore, the dimensionless wave
equation (57) becomes $u_{tt}=u_{xx}$ outside of that domain and up to its
front surface. Suppose that the experimental detector is positioned at $x=a$,
we want to propagate the data to a closer location $x=\tilde{a}$ with
$a<\tilde{a}<0$, where the zero point $x=0$ indicates the center of the
sandbox. The last endpoint of the sandbox is given at $x=\tilde{b}$ for
$\tilde{b}>0$. According to the absorbing boundary conditions mentioned in
Remark 3, the one-dimensional wave equation satisfies the left-traveling
waves, i.e. $u_{t}=u_{x}$ at $x=a$. Similarly, the boundary condition
$u_{t}=-u_{x}$ at $x=\tilde{b}$ is exactly satisfied by right-traveling waves.

###  5.1 Data preprocessing

In our previous works with the frequency domain data, we observed that the
signals’ quality was not good; cf. e.g. [25, Figure 3] and [14, Figure 3].
That challenge is easy to understand because of several unwanted factors
coming from, e.g., the surrounding items, the available Wi-Fi signal, etc. It
was ended up with the use of the data propagation that enables us to obtain
the near-field data with a much better quality. In the present study, we
experience the same and extend the applicability of the data propagation
technique to the case of time-dependent data. In particular, our question is
to find the boundary data $u\left(\tilde{a},t\right)$ for $a<\tilde{a}<0$ from
$u\left(a,t\right)=\varphi\left(t\right)$.

After the Fourier transform

$\hat{u}\left(x,k\right)=\mathop{\displaystyle\int}\limits_{0}^{\infty}u\left(x,t\right)e^{\text{i}kt}dt,$

the absorbing boundary condition becomes

$\hat{u}_{x}\left(x,k\right)+\text{i}k\hat{u}\left(x,k\right)=0,\quad
x\in\left[a,\tilde{a}\right].$

Therefore, one has

$\partial_{x}\left[\hat{u}\left(x,k\right)e^{\text{i}kx}\right]=0,$

which leads to

$\hat{u}\left(x,k\right)e^{\text{i}kx}=\hat{u}\left(a,k\right)e^{\text{i}ka}.$

Since we know the Dirichlet boundary condition measured at $x=a$, we find that

$\hat{u}\left(a,k\right)=\mathop{\displaystyle\int}\limits_{0}^{\infty}u\left(a,t\right)e^{\text{i}kt}dt=\hat{\varphi}\left(k\right).$

Thus, we obtain
$\hat{u}\left(x,k\right)=\hat{\varphi}\left(k\right)e^{\text{i}ka}e^{-\text{i}kx}$.
Then, using the inverse Fourier transform, we derive that for any
$x\in\left[a,\tilde{a}\right]$,

$u\left(x,t\right)=\mathop{\displaystyle\int}\limits_{0}^{\infty}\hat{\varphi}\left(k\right)e^{-\text{i}k\left(t+x-a\right)}dk=\varphi\left(t+x-a\right).$

Hence, $u\left(\tilde{a},t\right)=\varphi\left(t+\tilde{a}-a\right)$. This is
regarded as the time-shifted technique to move the data closer in a spatial
direction, given knowledge of the distance between the source position and the
front surface of the sandbox. Note that the operation $t+\tilde{a}-a$ is valid
because both time and space are in the dimensionless regime of (57).

The dimensions of the sandbox and the thickness of the styrofoam give us
$\tilde{b}=(0.35+0.07)/0.3=1.4$ and the propagated location
$\tilde{a}=-0.5/0.3\approx-1.67$, which is a bit far away from the sandbox’s
front surface covered with styrofoam. Cf. Figure 3, we can estimate the
furthest distance between the detector and the zero point centered in the
sandbox. Particularly, we compute

$\displaystyle a=$ $\displaystyle-\text{maximal
distance}\left(\text{antenna/sandbox }+\text{ front surface of sandbox/center
O }\right.$ $\displaystyle\left.+\text{
antenna/detector}\right)=-\left(1.7+\sqrt{0.42^{2}+0.25^{2}}+0.17+\sqrt{10^{2}+10^{2}}\right)/0.3$
$\displaystyle\approx$ $\displaystyle\;-8.33.$

Now we are in the position to complete the preprocessing procedure for our raw
far-field data.

  1. 1.

Background subtraction: For every location of the source, we subtract the
reference data from the far-field measured data. The reference data are those
we measure without the presence of the experimental object in the sandbox. It
is necessary to use this subtraction since the dielectric constant of dry sand
is around 4; cf. [10]. Currently, we measure the data twice since we do not
know how to deal with the signal distractions. The subtraction is proposed to
hope that it can reduce uncertain noises involved in our raw data. This is
also another technical difficulty we face when working on the experiment. A
similar procedure was implemented in [14, 18]. We hope to measure only once in
the future. But to do this, we need additional analytical and numerical
efforts.

  2. 2.

Data propagation: We apply the technique to $u(a,t)=\varphi(t)$ and obtain the
near-field data $u\left(\tilde{a},t\right)$ for each source location.

  3. 3.

Optimal detector: For each source position, our detector records signal in a
relatively small area. We found computationally that the optimal detector is
determined by the maximal value of $\left|u\left(\tilde{a},0\right)\right|$
with respect to all detector locations. Thus, for each source position, we use
only the data from this detector. Naturally, optimal detectors are different
for different positions of the source. Computationally, for each source
position, we note that $\left|u\left(\tilde{a},0\right)\right|\neq 0$ but
small, which is different from the theory that it should be 0, see Figure 4b.
The data after this step can be denoted by
$u_{\text{op}}\left(\tilde{a},t\right)$ for each source location.

  4. 4.

Data truncation: After choosing the optimal detector for each source position,
a truncation is applied to eliminate useless signals. For every time-dependent
signal, we can determine the location (in time) of the maximal absolute value
of the chosen signal. In other words, such a location is found by
$\max_{t\in[0,T]}\left\\{\left|u_{\text{op}}\left(\tilde{a},t\right)\right|\right\\}$.
Then all signals outside of the interval centered around that location (with
radius being 30 time samples) are deleted; see Figures 4a, 4b and 4c for a
typical sample of near-field data for our reference model. The data after this
step can be denoted by $u_{\text{op}}^{\text{tr}}\left(\tilde{a},t\right)$.

  5. 5.

The delay-and-sum procedure: The so-called delay-and-sum procedure is commonly
used in SAR imaging to preprocess the data [12]. We apply this procedure to
the absolute value of the truncated data
$u_{\text{op}}^{\text{tr}}\left(\tilde{a},t\right)$ and obtain
$u_{\text{sad}}\left(x_{s},y_{s}\right)$ for $x_{s}$ being the source variable
and $y_{s}$ being the location variable in the slant range plane. Then we
truncate those values which are less than 95% of the maximal value of the
delay-and-sum data. In other words, we obtain
$u_{\text{sad}}^{\text{tr}}(x_{s},y_{s})$ by

(58) $\displaystyle
u_{\text{sad}}^{\text{tr}}\left(x_{s},y_{s}\right)=\begin{cases}u_{\text{sad}}\left(x_{s},y_{s}\right)&\text{if
}u_{\text{sad}}\left(x_{s},y_{s}\right)\geq
0.95\max_{x_{s},y_{s}}\left(\left|u_{\text{sad}}\right|\right),\\\
0&\text{elsewhere}.\end{cases}$

The preprocessed data $u_{\text{sad}}^{\text{tr}}$ are used then in our
inversion procedure. See Figure 4d for the processed delay-and-sum data for
our reference model.

As to the data preprocessing mentioned above, we would like to remark that the
conventional SAR imaging software normally applies to two procedures to the
data $F(\mathbf{x_{0}},t)$ before getting $f_{0,n}$ for our inversion. That
procedure includes the use of the matched filtering and of the delay-and-sum
technique. However, our numerical experience in [21] was that after the delay-
and-sum alone, it provides a sufficiently good information about the location
on the slant range plane of the target represented by the slant range image.
This observation prompts one to just use the delay-and-sum technique and apply
the resulting data to solve the 1D CIP along the central line of the antenna.
Here, we do not detail the internal steps of the delay-and-sum procedure and
refer instead to section III-A of our first work [21]. We note, however, that,
compared with [21], a significantly new element of the delay-and-sum procedure
of this paper is that we apply it to absolute values of our (experimental)
preprocessed data, rather than to just real parts of the data in all previous
publications. The primary reason is that the experimental data are very
challenging, compared with the simulated data. By testing with the reference
model, Figure 5 is a particular showcase of why we should use the absolute
value of the data. In particular, we depict in Figures 5a and 5c the
illustrations of the delay-and-sum data when using, respectively, the real
value and the absolute value of the preprocessed raw data. Due to the
dimensions of the sandbox and the increment of the moving antenna (see Figures
2–3 and their captions), we know that the sandbox should be between the 10-th
and 20-th points of the source position. Then after the truncation (58), it is
clear that in terms of cleaning unnecessary data, the use of the absolute
value of the data works rather better than its real value; see Figures 5b and
5d. Henceforth, if one does not have a decent data preprocessing prior solving
the 1-D coefficient inverse problem, it is practically hard to obtain a good
reconstruction. The rest of formulas of the sum-and-delay procedure are the
same as in section III-A of [21].

(a)

(b)

(c)

(d)

Figure 4: A typical sample of the time-dependent experimental data for our
experimental model after the subtraction of the background data. (a) The
experimental data before data propagation. (b) The raw data after data
propagation. (c) The corresponding experimental data after truncation. (d) The
truncated delay-and-sum data. The delay-and-sum procedure was applied to the
absolute values of the data, which is a new element here.

(a)

(b)

(c)

(d)

Figure 5: (a) The delay-and-sum data when using the real part of the raw data.
(b) The truncated delay-and-sum data when using the real part of the raw data.
(c) The delay-and-sum data when using the absolute value of the raw data. (d)
The truncated delay-and-sum data when using the absolute value of the raw
data.

###  5.2 Some details

For our experimental result, we choose the calibration factor $CF=43.17$ and
parameters $\lambda=2.25,\beta=0.33,\gamma=10^{-10}$. Even though our
theoretical results are valid for large values of $\lambda$, it is observed
numerically that $\lambda\in[1,3]$ works well. The same interval for $\lambda$
was used in our past numerical studies; see e.g. [14, 16, 17, 19]. It is quite
natural that the $CF$ for the experimental data is different from the one for
computationally simulated data we did in [21].

(a)

(b)

Figure 6: (a) Photo of the experimental object. Its dimensions are about
$17\times 7\times 7$ in $\text{cm}^{3}$ (H$\times$W$\times$L). (b) The
computed image of (a), where the location is presented in the dimensionless
regime. The computed slant range dielectric constant is 23.9 and it was
established in [14] that the true dielectric constant of this object is 23.8.

We remark that $T$ is determined by two times of the maximal distance between
the line of sources $L_{\text{src}}$ and the furthest point of the sandbox. In
this sense, we take $T=20>2(\tilde{b}-a)\approx 19.46$ and consider 1000
points of time samples.

In Figure 4, we show how our data preprocessing works from the raw data to the
processed delay-and-sum data, considering the reference model. As can be seen
in Figure 4a, the data propagation is applied to avoid the direct signal from
the antenna, which is as strong as the one scattered by the target (see Figure
4b). Since we do not want any “parasitic” signal, the truncation, based on the
absolute value of the signal, is used, which leads to Figure 4c.

Note again that our significant novelty here is the use of the absolute value
of the data, rather than its real part in the conventional studies of SAR.
This use is because we observe numerically that it improves the location of
the slant range target in the slant range plane, see Figure 5. In the next
step, the delay-and-sum procedure with the truncation (58) is applied to the
resulting data, which finalizes our preprocessing procedure. Such data for the
inversion are depicted in Figure 4d.

Prior to comments on the numerical results, we would like to mention a
postprocessing procedure that we use after the inversion. From the truncated
delay-and-sum data, we are able to form a rectangular area that contains those
preprocessed data of $>95$% values of the maximal absolute value, see again
(58). We denote that rectangle by
$\left[l_{1},l_{2}\right]\times\left[s_{1},s_{2}\right]$, where $l_{j},s_{j}$
($j=1,2$) are standing for the location and source coordinates, respectively.
Tentatively, we denote the computed slant range function after the inversion
$\tilde{\varepsilon}_{\text{1}}(x_{s},y_{s})$. For each source position
$x_{s}$, we define $\tilde{\varepsilon}_{\text{2}}(x_{s},y_{s})$ as

$\tilde{\varepsilon}_{\text{2}}(x_{s},y_{s})=\begin{cases}\max_{y_{s}\in\left[l_{1},l_{2}\right]}\left(\tilde{\varepsilon}_{\text{1}}(x_{s},y_{s})\right)&\text{if
}\tilde{\varepsilon}_{\text{1}}<0.5\max_{y_{s}\in\left[l_{1},l_{2}\right]}\left(\tilde{\varepsilon}_{\text{1}(x_{s},y_{s})}\right),\\\
\tilde{\varepsilon}_{\text{1}(x_{s},y_{s})}&\text{elsewhere}.\end{cases}$

This step is to improve the dimensions of the reconstructed slant range image
since many truncation steps were applied in the data preprocessing. Then we
compute $\tilde{\varepsilon}_{\text{comp}}(x_{s},y_{s})$ as

$\tilde{\varepsilon}_{\text{comp}}(x_{s},y_{s})=\begin{cases}1&\text{outside
}\left[l_{1},l_{2}\right]\times\left[s_{1},s_{2}\right],\\\
\max_{\left[l_{1},l_{2}\right]\times\left[s_{1},s_{2}\right]}\left(\tilde{\varepsilon}_{2}(x_{s},y_{s})\right)&\text{inside
}\left[l_{1},l_{2}\right]\times\left[s_{1},s_{2}\right].\end{cases}$

We find that the objects’ signal is approximately at the 600-th point of the
time samples. We now consider the position of the objects in terms of the
distance between their central point and the central $O$ of the sandbox. Cf.
Figure 3b, the central $O$ of the sandbox is 35 cm away from its front
surface. As the data are propagated to $\tilde{a}=-1.67$, these time samples
range from $\tilde{a}$ to $\tilde{b}$, meaning that the distance between the
front and rear surfaces of the sandbox should be $\tilde{b}-\tilde{a}=3.07$.
Hence, the center of the object should be located
$0.6\times(\tilde{b}-\tilde{a})-\tilde{a}=0.172$ away from the central point
$O$ of the sandbox. Hence, multiplying this number by 0.3 to get back the
dimension regime, the center of the bottle is around 5 cm upwards of the
central $O$. Note that “upwards” means we are going through the central point
$O$ in the direction to the rear surface of the sandbox. The calculated
position is not so accurate for the experimental bottle of water, as it should
only be 21 cm downwards. We believe that the data truncation we have imposed
on the near-field data leads to this inaccuracy.

While the location in depth of our experimental object is somewhat acceptable,
we accentuate that its length is noticeably accurate. Since the buried object
is put horizontally, the length we mean here is actually the height in
reality. Counting the yellow part in Figure 6b in combination with the
knowledge that the increment of source positions is 5 cm, it reveals that the
length of the computed bottle is 35 cm, centered at the 15-th source position.
The true dimension should be 17 cm (cf. Figure 6a), which is smaller than the
computed one. However, we find it still useful in the detection of land mines
since having information about smaller sizes is useful in the battlefield.
Besides, the reconstruction is fully accurate that the center of the bottle
should be around the 15-th source position.

Cf. [14], the true dielectric constant of the watered bottle is
$\varepsilon_{\text{true}}=23.8$. Compared with this, we report that its
computed dielectric constant is $\tilde{\varepsilon}_{\text{comp}}=23.9$. At
the same time, when $\lambda=0$, we find that
$\tilde{\varepsilon}_{\text{comp}}=15.9$, which tells the importance of the
Carleman weight function presented in our reconstruction. Our numerical result
shows that we are capable of reconstructing experimentally buried targets,
whose dielectric constants are far away from the unity.

##  6 Summary

In our recent work [21], we have proposed, for the first time, a globally
convergent convexification numerical method for a fully nonlinear model of SAR
imaging. Some works on nonlinear inverse SAR problems used locally convergent
methods; cf. e.g. [13, 31]. In this work, we have presented more rigorous
results. From the theoretical standpoint, we have proved new Theorems 7 and 8,
which claim the global convergence of the gradient method of the minimization
of our convexified Tikhonov-like functional. This is stronger than our
previous results for the gradient projection method (cf. e.g. [5, 14, 17, 19,
20]) since obviously, the second one is harder to implement. Furthermore,
these theorems explain, for the first time, our consistent observation of many
previous works on the convexification that the gradient descent method works
rather than the gradient projection method.

In this second work, we have presented a new numerical result for a quite
challenging case of imaging of buried targets. This context has an obvious
application in standoff identification of antipersonnel land mines and
improvised explosive devices. The numerical result is for experimental data.
We have used only one target for a logistical reason since the data collection
process is quite time consuming. We have demonstrated that our proposed method
provides a quite accurate value of the dielectric constant of an explosive-
like target with an acceptable location and size.

## Acknowledgments

V. A. Khoa acknowledges Prof. Paul Sacks (Iowa, USA) for recent support of his
research career.

## References

  * [1] A. D. Agaltsov, T. Hohage, and R. G. Novikov, An iterative approach to monochromatic phaseless inverse scattering, Inverse Problems, 35 (2018), p. 024001, https://doi.org/10.1088/1361-6420/aaf097.
  * [2] F. Ahmad, M. G. Amin, and G. Mandapati, Autofocusing of through-the-wall radar imagery under unknown wall characteristics, IEEE Transactions on Image Processing, 16 (2007), pp. 1785–1795, https://doi.org/10.1109/tip.2007.899030.
  * [3] N. V. Alekseenko, V. A. Burov, and O. D. Rumyantseva, Solution of the three-dimensional acoustic inverse scattering problem. The modified Novikov algorithm, Acoustical Physics, 54 (2008), pp. 407–419, https://doi.org/10.1134/s1063771008030172.
  * [4] M. Amin, Through-the-wall Radar Imaging, CRC Press, Boca Raton, FL, 2011.
  * [5] A. B. Bakushinskii, M. V. Klibanov, and N. A. Koshev, Carleman weight functions for a globally convergent numerical method for ill-posed Cauchy problems for some quasilinear PDEs, Nonlinear Analysis: Real World Applications, 34 (2017), pp. 201–224, https://doi.org/10.1016/j.nonrwa.2016.08.008.
  * [6] L. Beilina and M. V. Klibanov, Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems, Springer US, 2012, https://doi.org/10.1007/978-1-4419-7805-9.
  * [7] A. L. Bukhgeim and M. V. Klibanov, Global uniqueness of a class of multidimensional inverse problems, Doklady Akademii Nauk SSSR, 260 (1981), pp. 269–272.  English translation: A. L. Buhgeim, M. V. Klibanov, Global uniqueness of a class of multidimensional inverse problems, Soviet Mathematics Doklady, 24, 244–247, 1981.
  * [8] S. A. Carn, Application of synthetic aperture radar (SAR) imagery to volcano mapping in the humid tropics: a case study in East Java, Indonesia, Bulletin of Volcanology, 61 (1999), pp. 92–105, https://doi.org/10.1007/s004450050265.
  * [9] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer New York, 2013, https://doi.org/10.1007/978-1-4614-4942-3.
  * [10] D. J. Daniels, Surface-penetrating radar, Institution of Electrical Engineers, London, 1996.
  * [11] D. J. Daniels, A review of GPR for landmine detection, Sensing and Imaging: An International Journal, 7 (2006), pp. 90–123, https://doi.org/10.1007/s11220-006-0024-5.
  * [12] M. Gilman, E. Smith, and S. Tsynkov, Transionospheric Synthetic Aperture Imaging, Springer International Publishing, 2017, https://doi.org/10.1007/978-3-319-52127-5.
  * [13] O. Karakus and A. Achim, On solving SAR imaging inverse problems using nonconvex regularization with a Cauchy-based penalty, IEEE Transactions on Geoscience and Remote Sensing, (2020), pp. 1–13, https://doi.org/10.1109/tgrs.2020.3011631.
  * [14] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. H. Nguyen, A. J. Sullivan, and V. N. Astratov, Convexification and experimental data for a 3D inverse scattering problem with the moving point source, Inverse Problems, 36 (2020), p. 085007, https://doi.org/10.1088/1361-6420/ab95aa.
  * [15] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. H. Nguyen, A. J. Sullivan, and V. N. Astratov, An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivity from experimental backscattering data, Inverse Problems in Science and Engineering, (2020), https://doi.org/10.1080/17415977.2020.1802447.  (published online of this journal).
  * [16] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen, Convexification for a three-dimensional inverse scattering problem with the moving point source, SIAM Journal on Imaging Sciences, 13 (2020), pp. 871–904, https://doi.org/10.1137/19m1303101.
  * [17] M. V. Klibanov and A. E. Kolesov, Convexification of a 3-D coefficient inverse scattering problem, Computers & Mathematics with Applications, 77 (2019), pp. 1681–1702, https://doi.org/10.1016/j.camwa.2018.03.016.
  * [18] M. V. Klibanov, A. E. Kolesov, and D.-L. Nguyen, Convexification method for an inverse scattering problem and its performance for experimental backscatter data for buried targets, SIAM Journal on Imaging Sciences, 12 (2019), pp. 576–603, https://doi.org/10.1137/18m1191658.
  * [19] M. V. Klibanov, A. E. Kolesov, A. Sullivan, and L. Nguyen, A new version of the convexification method for a 1D coefficient inverse problem with experimental data, Inverse Problems, 34 (2018), p. 115014, https://doi.org/10.1088/1361-6420/aadbc6.
  * [20] M. V. Klibanov, J. Li, and W. Zhang, Convexification for the inversion of a time dependent wave front in a heterogeneous medium, SIAM Journal on Applied Mathematics, 79 (2019), pp. 1722–1747, https://doi.org/10.1137/18m1236034.
  * [21] M. V. Klibanov, A. V. Smirnov, V. A. Khoa, A. J. Sullivan, and L. H. Nguyen, Through-the-wall nonlinear SAR imaging, IEEE Transactions on Geoscience and Remote Sensing, (2021), https://doi.org/10.1109/TGRS.2021.3055805.  (published online of this journal).
  * [22] M. V. Klibanov and N. T. Thành, Recovering dielectric constants of explosives via a globally strictly convex cost functional, SIAM Journal on Applied Mathematics, 75 (2015), pp. 518–537, https://doi.org/10.1137/140981198.
  * [23] A. E. Kolesov, M. V. Klibanov, L. H. Nguyen, D.-L. Nguyen, and N. T. Thành, Single measurement experimental data for an inverse medium problem inverted by a multi-frequency globally convergent numerical method, Applied Numerical Mathematics, 120 (2017), pp. 176–196, https://doi.org/10.1016/j.apnum.2017.05.007.
  * [24] A. Lechleiter and D.-L. Nguyen, A trigonometric Galerkin method for volume integral equations arising in TM grating scattering, Advances in Computational Mathematics, 40 (2013), pp. 1–25, https://doi.org/10.1007/s10444-013-9295-2.
  * [25] D.-L. Nguyen, M. V. Klibanov, L. H. Nguyen, and M. A. Fiddy, Imaging of buried objects from multi-frequency experimental data using a globally convergent inversion method, Journal of Inverse and Ill-posed Problems, 26 (2018), pp. 501–522, https://doi.org/10.1515/jiip-2017-0047.
  * [26] L. Nguyen, M. Ressler, and J. Sichina, Sensing through the wall imaging using the Army Research Lab ultra-wideband synchronous impulse reconstruction (UWB SIRE) radar, in Radar Sensor Technology XII, K. I. Ranney and A. W. Doerry, eds., SPIE, 2008, https://doi.org/10.1117/12.776869.
  * [27] R. Novikov, The $\bar{\partial}$-approach toapproximate inverse scattering at fixed energy in three dimensions, International Mathematics Research Papers, 2005 (2005), p. 287, https://doi.org/10.1155/imrp.2005.287.
  * [28] R. G. Novikov, Multidimensional inverse spectral problem for the equation $-\Delta\psi+(v(x)-{E}u(x))\psi=0$, Functional Analysis and Its Applications, 22 (1989), pp. 263–272, https://doi.org/10.1007/bf01077418.
  * [29] R. G. Novikov, An iterative approach to non-overdetermined inverse scattering at fixed energy, Sbornik: Mathematics, 206 (2015), pp. 120–134, https://doi.org/10.1070/sm2015v206n01abeh004449.
  * [30] O. of the Chief of Ordnance, Catalog of enemy ordnance materiel, no. N2228-E, World War II Operational Documents, 1945, http://cgsc.contentdm.oclc.org/cdm/ref/collection/p4013coll8/id/2758.
  * [31] B. Picard and E. Anterrieu, Comparison of regularized inversion methods in synthetic aperture imaging radiometry, IEEE Transactions on Geoscience and Remote Sensing, 43 (2005), pp. 218–224, https://doi.org/10.1109/tgrs.2004.841482.
  * [32] V. G. Romanov, Inverse Problems of Mathematical Physics, Walter de Gruyter GmbH & Co.KG, 2019.
  * [33] S. Rotheram and J. T. Macklin, Inverse Methods for Ocean Wave Imaging by SAR, in Inverse Methods in Electromagnetic Imaging, Springer Netherlands, 1985, pp. 907–930, https://doi.org/10.1007/978-94-009-5271-3_11.
  * [34] J. A. Scales, M. L. Smith, and T. L. Fischer, Global optimization methods for multimodal inverse problems, Journal of Computational Physics, 103 (1992), pp. 258–268, https://doi.org/10.1016/0021-9991(92)90400-s.
  * [35] G. A. Showman, Stripmap SAR, in Principles of Modern Radar: Advanced techniques, Institution of Engineering and Technology, 2012, pp. 259–335, https://doi.org/10.1049/sbra020e_ch7.
  * [36] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen, Convexification for a 1D hyperbolic coefficient inverse problem with single measurement data, Inverse Problems & Imaging, 14 (2020), pp. 913–938, https://doi.org/10.3934/ipi.2020042.
  * [37] A. V. Smirnov, M. V. Klibanov, A. J. Sullivan, and L. H. Nguyen, Convexification for an inverse problem for a 1D wave equation with experimental data, Inverse Problems, 36 (2020), p. 095008, https://doi.org/10.1088/1361-6420/abac9a.
  * [38] J. Sochacki, Absorbing boundary conditions for the elastic wave equations, Applied Mathematics and Computation, 28 (1988), pp. 1–14, https://doi.org/10.1016/0096-3003(88)90130-0.
  * [39] N. T. Thành, L. Beilina, M. V. Klibanov, and M. A. Fiddy, Imaging of buried objects from experimental backscattering time-dependent measurements using a globally convergent inverse algorithm, SIAM Journal on Imaging Sciences, 8 (2015), pp. 757–786, https://doi.org/10.1137/140972469.
  * [40] A. N. Tikhonov, Numerical Methods for the Solution of Ill-Posed Problems, Springer Netherlands, Dordrecht, 1995.