Numerical reconstruction for 3D nonlinear SAR imaging via a version of the convexification method

Synthetic-aperture radar (SAR) imaging is a commonly used technique in reconstructing images of surfaces of planets, and in detecting antipersonnel land mines and improvised explosive devices; cf. e.g. [2, 3, 7, 8] for some essential backgrounds concerning SAR imaging. The conventional SAR imaging is based on the linearization via the Born approximation [7]. While providing accurate images of shapes and locations of targets, the linearization is problematic to produce accurate values of their dielectric constants, see, e.g. [10]. However, since dielectric constants characterize constituent materials of targets, then it is obviously important to accurately compute them.

To address this problem, a new convexification based nonlinear SAR imaging, CONSAR, was proposed in [10] for the first time. Testing of CONSAR for both computationally simulated and experimentally collected data has shown that dielectric constants of targets are indeed accurately computed [10, 19]. However, these publications are only about obtaining 2D SAR images on the so-called “slant range” plane.  In this case both the point source and the radiating antenna move along an interval of a single straight line. Thus, unlike [10, 19], the current paper is about testing the CONSAR technique for 3D SAR images. By our experience, collection of experimental SAR data along just one line is time consuming [19]. Therefore, we consider here the case when the source and the detector concurrently move along only three (3) parallel lines. In principle, one should have many such lines of course to get accurate solutions. Hence, because of an obvious lack of the information content in the input data, it is unrealistic to expect that resulting SAR images would be accurate in terms of shapes of targets. This is why sizes of computed images of our targets are quite different from the true ones. Nevertheless, we accurately image dielectric constants and locations of targets. Only computationally simulated data are considered here.

In this paper, our target application is the non-invasive inspections of buried targets when one tries to image the buried or occluded target using backscattering data collected by e.g., flying devices or roving vehicles (cf. [2, 9]). Having knowledge of the dielectric constant of the target plays a vital role in the development of future classification algorithms that distinguish between explosives and clutter in military applications. Thus, it is anticipated that the knowledge of dielectric constants of targets should decrease the false alarm rate.

In fact, the conventional time dependent SAR data when the source and the detector move concurrently along an interval of a straight line, depend on two variables. On the other hand, we want to image 3D objects. Therefore, the data for the corresponding Coefficient Inverse Problem (CIP) are underdetermined. This is why we work here with the above mentioned more informative case of three parallel lines. Although these data are also underdetermined, we still hope to reconstruct some approximations of 3D images of targets. It is worth mentioning that solving underdetermined or non overdetermined CIPs is omnipresent in applications that we have mentioned above. This is one of the main challenges that researchers encounter in practice is finding a good approximation of the target’s dielectric constant. This is also a long standing research project that the corresponding author and his research team have been working on for almost a decade, see, e.g. [11, 12, 13, 14, 15, 16, 17] for mathematical models and numerical methods to identify the material of the buried targets.

Solving for the target’s dielectric constant in this scenario is a CIP for a hyperbolic PDE. It is well-known that CIPs are both severely ill-posed and nonlinear thus, requiring a careful designation of an appropriate numerical method. Very recently, we have proposed in [10, 19] a new convexification-based algorithm to construct a 2D image of the target’s dielectric constant using SAR data. In fact, the approach of [10, 19] means a nonlinear SAR imaging, as opposed to the conventional linear one. We point out that the computational approach of [10, 19] to the SAR data is a heuristic one. Nevertheless, it was shown in [10] that this approach works well numerically for computationally simulated data. It is more encouraging that it was demonstrated in [10, 19] that this approach works well for experimentally collected data. The same approach to the SAR data is used in this paper. Using multiple source locations running along an interval of a straight line, our computational approach is to solve a number of 1D CIPs.

To solve numerically each of our CIPs, we transform first the hyperbolic equation into a nonlocal nonlinear boundary value problem (BVP), based upon the original idea of [20]. An approximate solution of this BVP is then found by minimizing a suitable weighted Tikhonov-like cost functional, which is proved to be globally strictly convex. Its strict convexity is based on the presence of the Carleman Weight Function (CWF) in it. The CWF is the function which is used as the weight function in the Carleman estimate for the principal part of the differential operator involved in that BVP. Ultimately, the 2D image which we construct for the case of each of our three lines can be understood as the collection of 1D cross-sections of the so-called slant-range plane established by the “moving” source/antenna positions. To construct the 3D image, we combine those 2D ones, see sections 3 and 5 for details.

Conventional numerical methods for CIPs are based on minimizations of various least squares cost functionals, including the Tikhonov regularization functional, see, e.g. [4, 5, 6]. It is well known, however, that such a functional is typically non convex and has, therefore, multiple local minima and ravines, which, in turn plague the minimization procedure. The convexification method was initially proposed in [21, 22] with the goal to avoid that phenomenon of multiple local minima and ravines of conventional least squares cost functionals for CIPs. While publications [21, 22] are purely theoretical ones, the paper [23] has generated many follow up numerical studies of the convexification of this research group, see, e.g. the above cited publications of this research group and references cited therein. We also refer to the recently published book [24].

The convexification uses a significantly modified idea of the Bukhgeim–Klibanov method [25] of 1981. However, in [25] and many follow up publications, only proofs of global uniqueness and stability theorems for CIPs were studied, see, e.g. [26, 27, 28, 29, 30] and references cited therein. On the other hand, the convexification is going in the numerical direction: it is designed for constructions of globally convergent numerical methods for CIPs. This is unlike the conventional locally convergent numerical methods for CIPs, which are based on minimizations of those functionals.

As stated above, depending on the differential operator involved in the CIP, the convexification method is based on an appropriate Carleman Weight Function to construct a weighted Tikhonov-like cost functional. The Carleman Weight Function plays two roles. First, it convexifies the conventional Tikhonov functional. Second, it helps to control nonlinearities involved in the underlying CIP.

Definition 1. We call a numerical method for a CIP globally convergent if there is a theorem which claims that this method reaches at least one point in a sufficiently small neighborhood of the true solution of that CIP without any advanced knowledge of a small neighborhood of this solution.

The convexification is a globally convergent numerical method because one can prove the convergence of the scheme towards the true solution starting from any point located in a given bounded set of a suitable Hilbert space and the diameter of this set can be arbitrary large.

As to the moving sources, we remark that in the publication [16] a version of the convexification method was developed for the first time both analytically and numerically for a 3D inverse medium problem for the Helmholtz equation with the backscattering data generated by point sources running along a straight line and at a fixed frequency. Then this result was confirmed on experimental data in [17, 18]. In [16, 17, 18], we were working with the data in the frequency domain, whereas the SAR data are the time-dependent ones. In addition, in the case of the SAR data, both source and detector move concurrently along a straight line. On the other hand, in [16, 17, 18] the source moves along a part of a straight line, whereas the detector moves independently along a part of a plane. We also refer to works of Novikov and his group [31, 32] for a different approach to CIPs with the moving source.

We also add that our recent publications [10, 19] about CONSAR were focusing only on the 2-D reconstructions. Therefore, the main difference between these papers and the current one is that here we obtain a 3-D image of the target dielectric constant, assuming that the data are collected from multiple lines of sources.

Our paper is organized as follows. In sections 2, 3 and 4, we revisit our setting for the CIP imaging and our convexification-based method proposed in [19]. In section 5, we verify the 3-D numerical performance of our method using computationally simulated data. Finally, we present concluding remarks in section 6.

Prior to the statement of the forward problem, we discuss the setting of SAR imaging. We focus only on the stripmap SAR processing; see e.g. [33] for details of this setup. With the non-invasive inspection of buildings and land mines detection being in mind, we prepare to have a radiating antenna and a receiver in the SAR device. As the antenna radiates pulses of the time-resolved component of the electric wave field, the receiver collects the backscattering signal. Our SAR device moves along a straight line. In this work, we assume that it moves it along multiple parallel lines, which, we hope, may create a 3-D image of the target. It is worth mentioning that in our inverse setting in section 3 below, we assume that the antenna and the receiver coincide, and form a point, which moves along certain parallel lines. This assumption is a typical one in SAR imaging; cf. [7]. Indeed, in practice the antenna and the receiver are rather close to each other, and the distance between them is much less than the distance between the line over which they concurrently move and the targets of interest. However, when generating the SAR data in computational simulations, we solve a forward problem, in which we assume that the antenna is a disk.

We denote below $\mathbf{x}=(x,y,z)\in\mathbb{R}^{3}.$ For $z^{0}>0$, we consider the source location $\mathbf{x}^{0}=(x^{0},y_{j}^{0},z^{0})\in L_{j}^{\text{src}}$ for $j=\overline{1,j_{0}}$, $j_{0}\in\mathbb{N}^{\ast}$. By choosing a certain value of $j_{0}$ and varying $x^{0}$ in the interval $(-L,L)$ for $L>0$, we define a finite set of lines $L_{j}^{\text{src}}$ of sources/receivers as follows:

Therefore, each line of sources $L_{j}^{\text{src}}$ is parallel to the $x-$axis. The antenna and the receiver coincide and are located at the point $\mathbf{x}^{0}=\left(x^{0},y_{j}^{0},z^{0}\right)\in L_{j}^{\text{src}}.$ The number $y_{j}^{0}$ is the same for all points of $L_{j}^{\text{src}}$ and this number is changing as lines $L_{j}^{\text{src}}$ changes. The number $z^{0}$ is the same for all lines (2.1). Assuming that the parameter $x^{0}\in\left(-L,L\right)$ in (2.1), we actually assume that the antenna/receiver point characterized by $x^{0}$ moves along each line $L_{j}^{\text{src}},$ $j=\overline{1,j_{0}}.$

Let the number $R>0.$ Our domain of interest is

This domain consists of two parts. The lower part $\Omega\cap\left\{z<0\right\}$ mimics the sandbox containing the unknown buried object. Meanwhile, the upper part $\Omega\cap\left\{z>0\right\}$ is assumed to be the air part, where we move our SAR device. For simplicity, we assume that the interface between the air and sand parts is located at $\left\{z=0\right\}$. Let $\varepsilon_{r}(\mathbf{x})\in C^{1}\left(\mathbb{R}^{3}\right)$ be a function that represents the spatially distributed dielectric constant. We assume that

Physically, assumption (2.2) is reasonable. Cf. e.g. [17], we know that the dielectric constant of the air/vacuum is identically one. Meanwhile, for the other materials (for example, the dry sand which follows the non-invasive inspections of our interest), the dielectric constant is larger than the unity.

Let $T>0$ be the observation time. Practically, $T$ is found by doubling the distance from the line of sources to the furthest point of the sandbox. For every source location $\mathbf{x}^{0}\in L_{j}^{\text{src}}$ (see (2.1)), we consider the following forward problem:

where $n=n\left(\mathbf{x}\right),\mathbf{x}\in\partial\Omega$ is the unit outward looking normal vector on $\partial\Omega.$ In (2.3), $u(\mathbf{x},\mathbf{x}^{0},t)$ is the amplitude of a time-resolved component of the electric wave field. The function $\nu(\mathbf{x})=c_{0}/\sqrt{\varepsilon_{r}(\mathbf{x})}$ represents the speed of light in the medium with $c_{0}$ is the speed of light in the vacuum. Here, the source function $F$ is defined as

where $\text{i}=\sqrt{-1}$. This source function models well the linear modulated pulse, where $\omega_{0}$ is the carrier frequency, $\alpha_{0}$ is the chirp rate [7], and $\tau_{0}$ is its duration. The function $\chi_{D}$ in (2.4) represents the disk-shaped antenna of our interest with its radius $D>0$, thickness $\widetilde{D}>0$ and centered at a point $\mathbf{x}\in L_{j}^{\text{src}}$. We rotate the disk via a coordinate transformation $\mathbf{x}\mapsto\mathbf{x}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ defined by

where $\theta$ is the elevation angle. Hereby, the function $\chi_{D}$ is given by

We consider the so-called absorbing boundary conditions in (2.3). Mathematically, the solution to the acoustic wave equation is sought on the whole space. As introduced above, the computational domain $\Omega$ is, however, finite and the boundary condition should be equipped with the underlying PDE. Moreover, at this point, numerical simulations show the non-physical back reflections of the waves incident at the boundaries. Therefore, the absorbing boundary conditions are used to minimize these spurious reflections; cf. e.g. [34] and references cited therein.

Our application being in mind is to identify the material of the buried object. Therefore, our inverse problem aims to compute the spatial distribution of the dielectric constant $\varepsilon_{r}(\mathbf{x})$ involved in the coefficient $\nu(\mathbf{x})$ in the PDE of (2.3).

Consider the subspace $H_{0}^{k}(\mathcal{K})\subset H^{k}(\mathcal{K})$ for $k\in\mathbb{N}^{\ast}$. In this work, we introduce

As mentioned above, several transformations lead us to an overdetermined BVP. Since this problem is nonlinear, then any conventional optimization-based approach, including the Tikhonov regularization and the least-square method may suffer from the phenomenon of local minima and ravines. This is the main reason which has originally prompted the corresponding author to develop the convexification approach [21, 22]. The main ingredient of the method and its variants is using an appropriate Carleman Weight Function to “convexify” the cost functional. The associated high-order regularization terms usually play the role in controlling nonlinear terms involved in the PDE operator. Thereby, one obtains a unique minimizer and then, the global convergence of the gradient-like minimization procedure is attained. Depending on the principal parts of differential operators, one can choose different Carleman Weight Functions, see, e.g. [24] for a variety of Carleman estimates for partial differential operators. In this work, we consider the following Carleman Weight Function for the linear operator $\partial_{Y}^{2}-2\partial_{Y}\partial_{\tau}$ (see (3.18)):

Denoting

we come up with the following weighted Tikhonov-like functional $J_{\lambda,\gamma}:H^{4}(\mathcal{K})\rightarrow\mathbb{R}_{+}$:

We choose the high-order regularization term in $H^{4}(\mathcal{K})$ because of the embedding $H^{4}(\mathcal{K})\subset C^{2}(\overline{\mathcal{K}})$. This $C^{2}$ regularity is helpful in controlling our nonlinear term in (3.18); see [20, Theorem 2]. Following our convexification framework in, e.g., [16, 15, 20], we introduce the set $B\left(r,q_{0},q_{1}\right),$

in which $r>0$ is an arbitrary number. However, in our computational practice we use the $H^{2}\left(\mathcal{K}\right)-$norm, see section 5.1.

Minimization Problem. Minimize the cost functional $J_{\lambda,\gamma}\left(V\right)$ defined in (4.2) on the set $\overline{B\left(r,q_{0},q_{1}\right)}.$

We are now in the position to go over our analysis in [20, 36] of the convexification-based method for the nonlinear boundary value problem (3.18)–(3.21).

By the conventional regularization theory (cf. [37]), the existence of the ideal noiseless data $g_{0}^{\ast}\left(\tau\right),g_{1}^{\ast}\left(\tau\right)$ is assumed a priori. The existence of the corresponding exact coefficient $p^{\ast}\left(Y\right)$ and the exact function $w^{\ast}\left(Y,\tau\right)$ in (3.9)–(3.11) is assumed as well. Moreover, this function $p^{\ast}\left(Y\right)$ should also satisfy conditions $p^{\ast}(Y)=0$ for $Y<0$ and $Y>\sqrt{\overline{c}}$ as ones for $p(Y)$. Having the function $w^{\ast}\left(Y,\tau\right)$, we apply the above transformations (3.14), (3.17) to obtain the corresponding function $V^{\ast}\left(Y,\tau\right)$. Henceforth, it is natural to assume that functions $q_{0}^{\ast}$ and $q_{1}^{\ast}$ are the noiseless data $q_{0}$ and $q_{1}$ respectively; see (3.19) and (3.20). Thus, we assume below that the function $V^{\ast}\in B(r,q_{0}^{\ast},q_{1}^{\ast})$. Hence, it follows from (3.18)–(3.21) and (4.1) that

Consider a sufficiently small number $\sigma\in(0,1)$ that characterizes the noise level between the data $q_{0},q_{1}$ and $q_{0}^{\ast},q_{1}^{\ast}$. To obtain the zero boundary conditions at $\left\{x=0\right\}$ in (3.19) and (3.20), consider two functions $F\in B(r,q_{0},q_{1})$ and $F^{\ast}\in B(r,q_{0}^{\ast},q_{1}^{\ast})$. We assume that

Consider functions $W,W^{\ast},$

Besides, the function $W+F\in B(3r,q_{0},q_{1})$, $\forall W\in B_{0}(2r)$. We modify the cost functional $J_{\lambda,\gamma}\left(V\right)$ as

It follows from Theorem 2 that the functional $\widetilde{J}_{\lambda,\gamma}$ is also globally strict convex on the ball $\overline{B_{0}(2r)}$ for $\lambda\geq\lambda_{2}=\lambda_{1}\left(\mathcal{K},3r,\alpha\right)\geq\lambda_{1}\left(\mathcal{K},r,\alpha\right).$ Also, by Theorem 2, for each value of the parameter $\lambda\geq\lambda_{2},$ the functional $\widetilde{J}_{\lambda,\gamma}\left(W\right)$ has a unique minimizer $W_{\min,\lambda,\gamma}$ on the set $\overline{B_{0}(2r)}$ and

It follows from the proof of Theorem 5 of [20] that inequality (4.6) plays an important role in the proof of Theorem 3.

Now, we state a theorem about the global convergence of the gradient descent method of the minimization of the functional $J_{\lambda,\gamma}\left(V\right)$. Let $\omega\in(0,1)$ be the step size of this method. Fix an arbitrary number $\vartheta\in(0,1/3)$. We restrict the starting point, denoted by $V^{(0)}$, to be an arbitrary point in the set $B\left(\vartheta r,q_{0},q_{1}\right)\subset B(r,q_{0},q_{1})$. Then the gradient descent method reads as:

where $V^{(n)}$ approaches the minimizer $V_{\min,\lambda,\gamma}$ as $n$ is large. This scheme is well-defined because $J_{\lambda,\gamma}^{\prime}\left(V^{(n-1)}\right)\in H_{0}^{4}(\mathcal{K})$; see Theorem 4.2. Hence, functions $V^{(n)}$ satisfy the same boundary conditions as in (4.3) for all $n$.

Below, we assume that our minimizer $V_{\min,\lambda,\gamma}$ belongs to the set $B\left(\vartheta r,q_{0},q_{1}\right)$, which is an interior point of the set $B\left(r,q_{0},q_{1}\right)$. This enabled us in [19] to prove the global convergence of the gradient descent scheme (4.11). This assumption is plausible because by (4.9) the distance between the ideal solution $V^{\ast}$ and the minimizer $V_{\min,\lambda,\gamma}$ is small with respect to the noise level $\sigma\in(0,1),$ although only in the $H^{1}\left(\mathcal{K}_{\alpha,\mu,b}\right)-$norm rather than in the stronger norm $H^{4}\left(\mathcal{K}\right)$, see (4.3).

Combining Theorems 4.3 and 4.4, we can prove that the functions $V^{(n)}$ converge to the ideal solution $V^{\ast}$ as long as the noise level $\sigma$ in the data tends to zero. In addition, we obtain the convergence of the corresponding sequence $\left\{p^{(n)}\right\}_{n=0}^{\infty}$ towards the ideal function $p^{\ast}$. It works via the expression $p^{(n)}(Y)=4\partial_{Y}V_{Y}^{(n)}(Y,0)$ in (3.22).