**Preface to "Microwave Sensing and Imaging"**

Microwave sensing and imaging is acquiring an ever-growing importance in several applicative fields, such as non-destructive evaluations in industry and civil engineering, subsurface prospection, security, and biomedical imaging. Microwave techniques can be, in principle, used to retrieve information on some physical parameters of the inspected targets (dielectric properties, shape, etc.) by using safe electromagnetic radiations and cost-effective systems.

Although great technological advances have been attained in recent years, there is still a great deal of scientific research activity in this field, with the aim of further improving imaging systems and techniques. More efficient and reliable measurement systems are continuously being designed and validated on a case-by-case basis to address specific scenarios. Moreover, great attention has been paid to the development of effective data-processing algorithms, which are able to solve the underlying electromagnetic inverse scattering problem (which is generally nonlinear and ill-posed) in order to retrieve the required information on the inspected targets from the measured scattered-field samples. Finally, efficient forward solvers, which are fundamental for modeling the electromagnetic interactions between the interrogating fields and the targets, are proposed for both the development of the inversion approaches and for the definition/validation of the imaging setups.

This book, which reprints a Special Issue of the journal Sensors, provides some recent insights into microwave sensing and imaging systems and techniques, with reference to the topics outlined above.

> **Andrea Randazzo, Cristina Ponti, and Alessandro Fedeli** *Editors*

### *Article* **A Multi-Frequency Tomographic Inverse Scattering Using Beam Basis Functions**

**Ram Tuvi**

John A. and Katherine G. Jackson School of Geosciences, Institute for Geophysics, The University of Texas at Austin, Austin, TX 78758, USA; ram@ig.utexas.edu

**Abstract:** We present an overview of a beam-based approach to ultra-wide band (UWB) tomographic inverse scattering, where beam-waves are used for local data-processing and local imaging, as an alternative to the conventional plane-wave and Green's function approaches. Specifically, the method utilizes a phase–space set of iso-diffracting beam-waves that emerge from a discrete set of points and directions in the source domain. It is shown that with a proper choice of parameters, this set constitutes *a frame* (an overcomplete generalization of a basis), termed "beam frame", over the entire propagation domain. An important feature of these beam frames is that they need to be calculated once and then used for all frequencies, hence the method can be implemented either in the multifrequency domain (FD), or directly in the time domain (TD). The algorithm consists of two phases: in the processing phase, the scattering data is transformed to the beam domain using windowed phase–space transformations, while in the imaging phase, the beams are backpropagated to the target domain to form the image. The beam-domain data is not only localized and compressed, but it is also physically related to the local Radon transform (RT) of the scatterer via a local Snell's reflection of the beam-waves. This expresses the imaging as an inverse local RT that can be applied to any local domain of interest (DoI). In previous publications, the emphasis has been set on TD data processing using a special class of localized space–time beam-waves (wave-packets). The goal of the present paper is to present the imaging scheme in the UWB FD, utilizing simpler Fourier-based data-processing tools in the space and time domains.

**Keywords:** inverse scattering; imaging; wave propagation; beam summation methods

#### **1. Introduction**

Inverse scattering deals with determining the shape and the composition of an unknown object from measurements of the scattering field data due to a known illumination. This area has a wide range of medical, geophysical, oceanographical, industrial, etc., applications, using electromagnetic, acoustic, elastic, or seismic waves [1–5]. Inverse scattering problems are, in general, non-linear and highly ill-posed, hence accurate solutions typically require iterative schemes and are limited to rather small configurations in the order of wavelengths. For large domains, practical algorithms rely on linearized weak scattering formulations using the Born, Rytov, Physical optics, or other single scattering approximations [2,5] which linearize the relation between the target and the field and provide the basis for diffraction tomography (DT) reconstruction [6].

Inverse scattering requires diversity and relies on the wave data in hand. Depending on the application, it may involve multiple excitation frequencies (or short-pulse response) and/or several illumination directions. With the overall complexity of the problem, full utilization of the wave data is essential to formulate an efficient, robust, and accurate algorithm.

Beam summation (BS) methods are when the wave field is expressed as a superposition of collimated beam waves. Here, we use the generic term "beam" for both the FD formulations where the propagators are Gaussian beams, and for the TD formulations

**Citation:** Tuvi, R. A Multi-Frequency Tomographic Inverse Scattering Using Beam Basis Functions. *Sensors* **2022**, *22*, 1684. https://doi.org/ 10.3390/s22041684

Academic Editors: Andrea Randazzo, Cristina Ponti and Alessandro Fedeli

Received: 2 December 2021 Accepted: 3 February 2022 Published: 21 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where the propagators are pulsed beams. This provides the proper physical basis for such robust DT reconstruction. Like the plane-wave (PW) spectrum approach, the BS provides a complete spectral basis which is required for DT reconstruction, yet unlike plane-waves, the beam waves provide spatial resolution and they can easily and efficiently be propagated in inhomogeneous media along ray trajectories. Unlike rays, on the other hand, they provide a uniform spectral basis and are insensitive to geometrical optics transition regions. Thus, beams provide a way to convert the wave problem to a ray-based skeleton.

BS methods can be classified into two classes (see review in [7,8]). The first class addresses radiation from localized (point) sources by expressing the field as an angular superposition of collimated beams that emerge radially from the source. This formulation has been derived asymptotically [9,10], but later on it was formulated as an exact spectral identity using complex source beams [11,12] and was also extended to the time domain (TD) [12]. Consequently, it has been used in various applications of propagation, scattering, and inverse scattering.

The other class addresses radiation from extended (aperture) sources, where the field is expressed as a sum of beam propagators that emerge from a discrete phase–space lattice of points and directions in the aperture. These formulations utilize local window (e.g., Gaussian) functions to transform the data to the beam domain, and then propagate the data using beams. These formulations are utilized in the present work where we analyze the data on the measurement aperture and then backpropagate it to the scatterer domain. Early implementation of this approach was based on a Gabor series expansion of the planar field [13–16] and therefore suffered from two major drawbacks: (a) The Gabor expansion coefficients (the beam amplitudes) are notoriously non-local and unstable, and (b) the beam lattice is frequency-dependent, hence a new lattice should be calculated at each frequency. Both difficulties have been mitigated in the ultra-wide band phase–space beam summation (UWB-PS-BS) method [17,18] which utilizes the overcomplete windowed Fourier transform (WFT) frames. In linear algebra, a frame of an inner product space is a generalization of a basis of a vector space to overcomplete sets. In signal processing, frames provide a redundant, stable way of representing a signal, instead of the Gabor series. The formulation is structured upon a *frequency independent* lattice of beams that emerge from a discrete phase–space lattice of points and directions in the aperture, and utilizes iso-diffracting Gaussian beams (ID-GB) whose propagation parameters are *frequency independent.* These properties make this representation efficient for wideband applications and also allow an extension of the theory to the time-domain (TD) [19].

A major step forward has been the proof in [20] that these phase–space sets of beams constitute a frame not only in the aperture domains but actually everywhere in the propagation domain. This implies that these beam basis functions may be used to expand not only the sources and the field, but also any function of space and in particular the medium inhomogeneities, a property that is being used in our beam-based inverse scattering theory. The theory has been proved originally in the frequency domain (FD) [20] and then extended to the TD in [21].

As noted above, it has been recognized long ago that BS provides the proper physical basis for a robust DT reconstruction. Examples for point-sources configurations can be found in [22–30]. For configuration where the data is measured over a wide aperture (see Figure 1a), it is more suitable to use the extended sources approach noted above: see [31–37] and [38–44] for medium reconstruction over a homogenous or inhomogeneous background, respectively. Without going into detailed comparison, the main difference between the methods is in the data representation phase (see [42] for a detailed comparison).

**Figure 1.** Diffraction tomography and the *K*-space diffraction tomography identity in the spatial and spectral domains. (**a**) The physical configuration. The unknown object *O*(**r**) is located between two measurement planes. At *z* = *z*<sup>1</sup> we have an array of sources/receivers, while at *z* = *z*<sup>2</sup> we have an array of receivers. The object is illuminated by the plane-wave (black arrows) of Equation (3) propagating in the direction ◦ *κi* . In red we plot pulsed plane-wave illumination of Equation (3). The scattered field is measured on the *zj* planes. (**b**) The DT identity of (8). The plane-wave spectrum of the scattered field ˆ*u*˜*<sup>s</sup> j* (*ξ*) is mapped to values of *O*¯(**K**) over a shifted Ewald sphere **K** = *k*( ◦ *<sup>κ</sup><sup>j</sup>* <sup>−</sup> ◦ *κi* ). The data measured on the *j* = 1, 2 planes are illustrated by red dashed and blue solid-line hemispheres, respectively.

In this work we review the beam-based local inverse scattering theory derived in [36,37], which is based on the frame-based UWB-PS-BS theory discussed above. As noted there, the theory is structured on a frequency-independent phase–space sets of beams that constitute frames everywhere in the propagation domain. This beam frame formulation enables the expansion of both the medium inhomogeneities and the scattering data with the same set of beam-basis functions, thus enabling a direct inversion over the beam domain. In previous publications, the emphasis has been set on TD data processing using special localized space–time beam-waves (wave-packets). This requires somewhat sophisticated mathematics and processing tools. In the present paper, on the other hand, we utilize a simpler FD Fourier-based data-processing approach followed by an integration over the relevant frequency band. The paper makes extensive references to equations or figures in [36,37]. Therefore, to simplify the presentation, we refer to them by the prefixes I and II, respectively.

The advantages of our beam-based DT over the conventional DT approach are:


*e*. The beam-based imaging enables backpropagation and imaging over a non-homogeneous background.

The presentation below starts in Section 2 with a review of the main concepts in DT. We proceed in Section 3 by reviewing elements of the beam representation, and in particular of the UWB-PS-BS and the BF concept. The beam-based DT is presented in Section 4, where, as discussed above, we emphasize the multi-frequency data processing as opposed to the more complicated TD processing used in [36,37]. The presentation ends in Section 5 with a practical description of the algorithm, including the choice of the various parameters and numerical examples.

#### **2. UWB Diffraction Tomography in the Spectral Plane-Wave Domain: A Review**

This section reviews the conventional plane-wave DT algorithms. Referring to the configuration in Figure 1a, we consider the two alternative schemes: The *angular diversity* scheme where the data is measured for several illumination directions ◦ *κ<sup>i</sup>* at a given frequency (Section 2.3), and the *frequency diversity* scheme where the data is measured over a wide frequency band for a single illumination direction ◦ *κ<sup>i</sup>* (Section 2.4). The frequency diversity scheme may also be performed using a short-pulse illumination and calculated directly in the TD [45] (see also I– Section 3-B in [36]).

#### *2.1. Problem Description—Physical Configuration*

The physical configuration is illustrated in Figure 1a, where the object is located between two measurement planes, at *z* = *z*<sup>1</sup> < 0 and at *z* = *z*<sup>2</sup> > 0. We assume a 3*D* coordinate frame **r** = (**x**, *z*) where the z-coordinate is normal to the measurement planes, and **x** = (*x*1; *x*2) are the transversal coordinates. The data is collected over a wide frequency band Ω ∈ [*ω*min, *ω*max]. The theory is presented here in the FD, but we also discuss the TD formulation for completeness and clearer interpretation. Field constituents in these domains are related via the temporal Fourier transform

$$
\hat{u}(\omega) = \int dt \, u(t) \, e^{i\omega t} \,\tag{1}
$$

where FD constituents are tagged by an over-hatˆ.

The unknown object is embedded in a uniform background wavespeed *v*<sup>0</sup> and assumed to be lossless and non-dispersive. It is described by the unknown wavespeed *v*(**r**), and we define the "object function"

$$O(\mathbf{r}) = \left(v\_{\mathbb{0}}/v(\mathbf{r})\right)^2 - 1\tag{2}$$

(see Equation (7)).

The scattering data may be collected as a function of frequency using time-harmonic plane-wave excitation, or directly in the TD utilizing short-pulse plane-wave. These excitations are given by

$$\mathfrak{u}^{i}(\mathbf{r},\omega) = \mathbb{P}(\omega)e^{i\mathbf{k}\check{\mathbf{x}}^{i}\cdot\mathbf{r}}, \quad \mathfrak{u}^{i}(\mathbf{r},t) = F(t - \upsilon\_{0}^{-1}\check{\mathbf{x}}^{i}\cdot\mathbf{r}), \tag{3}$$

where *F*ˆ(*ω*) is the source spectrum and *k* = *ω*/*v*<sup>0</sup> is the wavenumber. The incident wave propagates in the direction

$$\hat{\mathbf{x}}^{\hat{i}} = (\mathfrak{f}^{\dot{i}}, \mathbb{Q}^{\dot{i}}) = \sin \theta^{\dot{i}} \cos \phi^{\dot{i}} \hat{\mathbf{x}}\_1 + \sin \theta^{\dot{i}} \sin \phi^{\dot{i}} \hat{\mathbf{x}}\_2 + \cos \theta^{\dot{i}} \hat{\mathbf{z}},\tag{4}$$

with (*θ<sup>i</sup>* , *φ<sup>i</sup>* ) being the polar angles with respect to the *z* axis, and over-circles denote unit vectors.

The scattered fields measured over the *zj* planes, *j* = 1, 2 are denoted as *u*ˆ*<sup>s</sup> j* (**x**, *ω*) (see Figure 1a). The PW spectral representation of *u*ˆ*<sup>s</sup> j* (**r**) is defined via

$$\hat{u}\_{\vec{j}}^{s}(\mathfrak{F}) = \varepsilon^{\pm ik\tilde{k}\cdot z\_{\vec{j}}} \int d^2 \mathbf{x} \,\hat{u}\_{\vec{j}}^{s}(\mathbf{x}) e^{-ik\tilde{\xi}\cdot \mathbf{x}}, \quad \tilde{\mathfrak{J}} = \sqrt{1 - \mathfrak{F}\cdot \mathfrak{F}}, \quad \text{Im}\,\tilde{\mathfrak{J}} \ge 0. \tag{5}$$

where lower and upper signs correspond to *j* = 1 and 2, respectfully. We added the *e* ±*ikζzj* phase term in (5) in order to normalize the spectral PWs to the *z* = 0 plane instead of *z* = *zj* planes.

Note that we use here the frequency normalized spectral coordinates *ξ* = **k***x*/*k* which are related to the PW direction via *ξ* = (*ξ*1, *ξ*2) = sin *θ*(cos *φ*, sin *φ*), where (*θ*, *φ*) are the conventional spherical angles with respect to the *z*-axis so that the scattered PWs propagate in the unit vector direction

$$\stackrel{\circ}{\kappa}\_{\circ} = (\not\xi, \mp\zeta) = (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta), \quad \not\jmath = 1, 2. \tag{6}$$

The spectral ranges |*ξ*| < 1 and |*ξ*| > 1 define the *propagation spectrum* and *evanescent spectrum*, respectively. Typically, DT formulations are restricted only to the propagation spectrum data (see discussion after Equation (9)).

#### *2.2. The DT Identity*

According to the weak scattering (first Born) approximation of the Lippmann–Schwinger integral equation, the scattered field can be expressed as [5]

$$\hat{u}^s(\mathbf{r}) = k^2 \int\_V d^3r' \hat{u}^i(\mathbf{r'}) \mathcal{O}(\mathbf{r'}) \mathcal{G}(\mathbf{r}, \mathbf{r'}),\tag{7}$$

where *<sup>G</sup>*<sup>ˆ</sup> = *<sup>e</sup>ik*|**r**−**r**| <sup>4</sup>*π*|**r**−**r**| is the 3D Green's function in the uniform background. This approximation is valid if *O*(**r**) 1 and in addition *kLO*max < 1, where *L* is the spatial support of *O* and *O*max is its maximal value.

Inserting (7) into (5) and using the spectral representation of *G*ˆ, we obtain (I-7),

$$\hat{\mu}\_{\rangle}^{s}(\boldsymbol{\xi}) \simeq \frac{k}{-2i\zeta} \bar{O}(\mathbf{K}) \Big|\_{\mathbf{K} = k(\stackrel{\circ}{\kappa}\_{\rangle} - \stackrel{\circ}{\kappa}\_{\}} \qquad |\boldsymbol{\xi}| < 1,\tag{8}$$

where ◦ *κ<sup>j</sup>* and ◦ *κ<sup>i</sup>* are given by (6) and (4), and

$$O(\mathbf{K}) = \int d^3r \, O(\mathbf{r}) e^{-i\mathbf{K}\cdot\mathbf{r}}, \qquad \mathbf{K} = (K\_1, K\_2, K\_z), \tag{9}$$

is the *K*-space distribution of *O*(**r**). Equation (8) is referred to as the *DT identity*. It relates the scattering data in the ◦ *κ<sup>j</sup>* directions to values of *O*¯(**K**) at the points **K** = *k*( ◦ *<sup>κ</sup><sup>j</sup>* <sup>−</sup> ◦ *κi* ). As illustrated in Figure 1b, these points define a *K*-space sphere with radius *k* that is centered at **K** = −*k* ◦ *κi* , which is referred to as the *shifted Ewald sphere*. Note from (6) that the left and right hemispheres (plotted as red and blue, respectively) correspond to data from the *zj* measurement plane with *j* = 1, 2, respectively.

The DT identity above applies only to the propagation spectrum |*ξ*| < 1. Adding the evanescent spectrum may improve the resolution. However, the contribution of the evanescent spectrum is exponentially weak and hence has a low signal to noise ratio. In addition, backpropagating this data to form the image amplifies the noise level exponentially. For these reasons, the evanescent spectrum contribution is usually neglected except for near field imaging schemes.

In view of the DT identity, one may obtain a full *K*−space coverage of the object function by measuring the scattering response for several illumination directions or several frequencies [2,5]. These alternative schemes are reviewed in the following sections.

#### *2.3. Object Reconstruction via Angular Diversity*

The angular diversity approach is illustrated in Figure 2a. Changing the illumination directions ◦ *κ<sup>i</sup>* while keeping the operational frequency *k* constant changes the centers of the shifted Ewald sphere and provides a different coverage of the *K* space. Aggregating the response for several illumination directions recovers *O*¯(**K**). Note that for lossless (real) objects, *<sup>O</sup>*¯(**K**) = *<sup>O</sup>*¯ <sup>∗</sup>(−**K**), so that only half of the K-space needs to be recovered.

**Figure 2.** *K*-space reconstruction. (**a**) Angular diversity reconstruction: Changing the direction of illumination ◦ *<sup>κ</sup><sup>i</sup>* and measuring the transmitted data only provides a coverage of the *<sup>K</sup>*−space within a sphere of radius <sup>√</sup>2*<sup>k</sup>* [5]. (**b**) Frequency diversity reconstruction: Changing the excitation frequency for a single illumination direction ◦ *<sup>κ</sup><sup>i</sup>* provides coverage of the *<sup>K</sup>*<sup>−</sup> space as indicated. The figure is plotted for ◦ *κ<sup>i</sup>* = ◦ *z*.

As one observes from Figure 2a, the transmitted data on *z*<sup>2</sup> recovers the *K*-space distribution of the object in a sphere of radius <sup>√</sup>2*<sup>k</sup>* about the origin. Thus, the object may be recovered using only transmitted data, as long as *k* is chosen to be large enough to provide full coverage of *<sup>O</sup>*¯(**K**). Note that in the limit of *<sup>k</sup>* <sup>→</sup> <sup>∞</sup>, the transmitted data hemispheres in Figure 2a reduce to planar surfaces normal to ◦ *κ<sup>i</sup>* that pass through the origin, thus providing the *K*-space representation of conventional X-ray tomography [46].

One option to reconstruct *O*(**r**) is to recover *O*¯(**K**) and then apply the inverse Fourier transform of (9). This approach requires interpolation of the data from the shifted Ewald spheres to a Cartesian *K*-domain grid [47], and therefore requires a large number of illuminations.

The "filtered backpropagation" reconstruction algorithm [6,47] overcomes this difficulty by circumventing the need to recover *O*¯(**K**) and operating, instead, directly on the scattering data. In this approach, the scattering data is multiplied by a spectral filter and then back-propagated to the object domain. This reconstruction approach is analogous to the X-ray tomography filtered backprojection algorithm of [48], where the filtered data is back-projected along straight lines.

#### *2.4. Object Reconstruction via Frequency Diversity (UWB Tomography)*

In the frequency diversity approach, the data is collected over a wide frequency spectrum Ω ∈ [*ω*min, *ω*max] for a fixed illumination direction. This can be done either in a frequency by frequency approach or by using a short-pulse illumination. As noted earlier, in this paper we emphasize the multi-frequency approach. The readers are referred to I–Section 3-B in [36] for the TD formulation, which has an important cogent physical interpretation, but is not utilized here.

As illustrated in Figure 2b for illumination along the positive *z* axis, changing the illumination frequency changes the radius of the shifted Ewald sphere. One observes that the reflected data recovers the *K*-space distribution of *O* inside a *π*/4 cone with an axis along the negative *Kz* axis and base radii between *k*min and *k*max, while the transmitted data recovers the complementary *K*-space part. As noted before, only half of the *K*-space is needed to recover the real function *O*. Otherwise, several illumination directions are required. More illuminations also add robustness.

As follows from Figure 2b, the reflection data on the *z*<sup>1</sup> plane mainly recovers the object variations along the *z* axis, while the transmitted data on the *z*<sup>2</sup> plane recovers the transversal variation (see also Snell's law interpretation in I-Section 3-B in [36]. Thus, for quasi-stratified media with weak transversal variations, it may be sufficient to measure only the reflected data on the *z*<sup>1</sup> plane, but may not be sufficient for objects with a large transversal variations. Another limitation is the missing data for |*Kz*| < 2*k*min (see Figure 2b), while the missing data for |*Kz*| > 2*k*max can be measured by using higher frequencies. As follows from Figure 2b, several illumination directions may increase the transversal resolution and also add data at small |**K**|. Note also that the maximal axial resolution for the case of normal incidence is *δ<sup>z</sup>* = *π*/*k*max.

The object can be reconstructed using an inverse transform of *O*¯(**K**). However, for the same reasons discussed in Section 2.3, a filtered backpropagation approach is preferable. Backpropagation can be calculated in several alternative ways. For simplicity, we present here the spectral integration approach. Given the scattering data *u*ˆ*s*(**x**, *ω*) over the *zj* planes, the backpropagated fields into the *z* > *z*<sup>1</sup> and *z* < *z*<sup>2</sup> regions are given by (see (5))

$$\text{ch}\_{\mathfrak{f}}^{b}(\mathbf{r}) = \left(\frac{k}{2\pi}\right)^{2} \int\_{|\mathfrak{f}|<1} d^{2}\xi^{\mathfrak{a}} \hat{\mathfrak{a}}\_{\mathfrak{f}}^{s}(\mathfrak{f}) e^{i\mathbf{k}\cdot(\mathfrak{f}\cdot\mathbf{x}\mp\zeta z)},\tag{10}$$

where we restrict the integration to the visible spectrum |*ξ*| < 1.

The "imaging field," or the "filtered backpropagated field" corresponding to the data on the *j* = 1, 2 plane is given by (II-2)

$$
\hat{I}\_{\rangle}(\mathbf{r},\omega) = \upsilon\_0^{-1} k^{-2} \hat{\mathbf{x}}^{\hat{i}} \cdot \nabla \left[ e^{-ik\hat{\mathbf{x}}^{\hat{i}} \cdot \mathbf{r}} \hat{u}\_{\rangle}^{b}(\mathbf{r},\omega) \right]. \tag{11}
$$

The corresponding "partial images" are obtained by summing over the relevant frequency band (II-3)

$$\mathcal{O}(\mathbf{r}) = 2 \text{Re} \, \frac{1}{\pi} \int\_0^\infty d\omega \, \mathbf{\hat{l}}\_{\mathbf{\hat{j}}}(\mathbf{r}, \omega). \tag{12}$$

If the data is given on both planes, then the "complete image" is given by

$$
\bullet \bullet (\mathbf{r}) = \bullet\_1(\mathbf{r}) + \bullet\_2(\mathbf{r}).\tag{13}
$$

The features of *O*(**r**) that are described by *O*˘ *<sup>j</sup>* have been discussed above in connection with Figure 2b. As noted there, in many situations it is sufficient to recover only *O*˘ 1.

The derivation of the filtered backpropagation imaging algorithm in (10)–(12) is done by inserting the Born approximated data of (8) into (10).

#### **3. Mathematical Background on the UWB-PS-BS**

#### *3.1. The Windowed Fourier Transform (WFT) Frame Representation of the Field*

As noted in the introduction above, the phase–space beam summation representation is based on the theory of WFT frame expansion of the field. Following [17], the theory is presented here in the context of radiation into the half-space *z* > 0 in a 3D coordinate space **r** = (**x**, *z*), **x** = (*x*1, *x*2), due to a time harmonic field *u*ˆ0(**x**, *ω*) defined over the plane *z* = 0. The WFT frame set {*ψ*ˆ*μ*(**x**, *<sup>ω</sup>*)} is defined by (Equation (22)] [17])

$$
\hat{\psi}\_{\mu}(\mathbf{x}, \omega) = \hat{\psi}(\mathbf{x} - \mathbf{m}\bar{\mathbf{x}})e^{i\mathbf{k}\mathbf{n}\_{\sharp}^{\mathrm{T}} \cdot (\mathbf{x} - \mathbf{m}\bar{\mathbf{x}})},\tag{14}
$$

with *ψ*ˆ being a localized window function (typically a Gaussian, see more details below) and *μ* = **m**, **n** being a 4-index. The frame elements are localized about the spatial (**x**) and spectral (*ξ*) phase space lattice

$$(\mathbf{x\_{m}}, \mathbf{\tilde{y}\_{n}}) = (\mathbf{m}\mathbf{\bar{x}}, \mathbf{n}\mathbf{\tilde{y}}) = (m\_1 \mathbf{\bar{x}}, m\_2 \mathbf{\bar{x}}; m\_1 \mathbf{\bar{y}}, m\_2 \mathbf{\bar{y}}),\tag{15}$$

where (*x*¯, ¯ *ξ*) defines the lattice unit cell. As will be shown, the points (**xm**, *ξ***n**) define the beams' initiation points and propagation directions (see Equation (21) below). To constitute a frame, the set above needs to fully cover the phase space, i.e., the unit cell area should be less than 2*π*, implying that

$$k\ddot{\mathbf{x}}\ddot{\mathbf{j}} = 2\pi\nu\_{\prime} \tag{16}$$

with *ν* < 1 being the overcompleteness parameter and the limit *ν* = 1 define the critically complete limit. As will be shown in Section 3.2 below, the frame over completeness provides a local and stable representation of the field (Equation (13) [17]), with it also being used to derive an UWB representation of the field (see Equations (Equations (33)–(35) [17]))).

The WFT frame can be used to expand *u*ˆ0(**x**) in the form

$$\mathfrak{A}\_0(\mathbf{x}) = \sum\_{\mu} \mathfrak{a}\_{\mu} \,\hat{\psi}\_{\mu}(\mathbf{x}). \tag{17}$$

In view of the overcompleteness, the coefficients set *a*ˆ*μ* is not unique. A particularly convenient set with a minimum -<sup>2</sup> norm is obtained by using the dual frame *ϕ*ˆ*μ*(**x**) which has the same structure as *ψ*ˆ*μ* in (15) except that the mother window *ψ*ˆ(**x**) is replaced by the dual mother window *ϕ*ˆ(**x**). The resulting canonical coefficient set is given by (Equation (23) [17]).

$$\mathfrak{A}\_{\mu} = \langle \mathfrak{A}\_{\mathbb{O}}(\mathbf{x}), \mathfrak{\phi}\_{\mu}(\mathbf{x}) \rangle \tag{18}$$

where *f* , *g*  = *f g*<sup>∗</sup> is the conventional L<sup>2</sup> inner product in the transverse coordinate **x**. The canonical coefficients *a*ˆ*<sup>μ</sup>* in (18) are readily identified as the local spectrum of *u*ˆ0(**x**) windowed with respect to *ϕ*ˆ*<sup>μ</sup>* about the phase–space points **xm**, *ξ***<sup>n</sup>** .

Generally, *ϕ*ˆ should be calculated numerically, for a given *ψ*ˆ and lattice *x*¯, ¯ *ξ* . However, if the lattice is sufficiently overcomplete, (*ν* - 1/3) *ϕ*ˆ ∝ *ψ*ˆ can be approximated by (Equation (11) [17])

$$
\phi(\mathbf{x}) \approx \nu^2 \hat{\psi}(\mathbf{x}) / \|\psi\|^2. \tag{19}
$$

There are mainly two reasons to prefer the use of this highly overcomplete parameter regime, even though it implies a larger number of terms in the phase–space expansion (17): (*i*)—as follows from (19), in this parameter regime *ϕ*ˆ is localized both spatially and spectrally, so that the expansion (17) comprises local and stable coefficients. (*ii*)—*ϕ*ˆ is given analytically via (19) and does not have to be to calculated numerically. Reason (*ii*) is critical for UWB problems where *ϕ*ˆ needs to be found for each *ω*.

The radiated field in *z* > 0 is obtained now by replacing *ψ*ˆ*μ*(**x**) in Equation (17) by beam propagators (Equation (24) [17])

$$\mathfrak{A}(\mathbf{r}) = \sum\_{\mu} \mathfrak{A}\_{\mu} \hat{\Psi}\_{\mu}^{+}(\mathbf{r}),\tag{20}$$

$$\hat{\Psi}^+\_{\mu}(\mathbf{r}) = \left(\frac{k}{2\pi}\right)^2 \int d^2\xi \,\hat{\Psi}\mu(\xi)e^{ik(\xi\cdot\mathbf{x}+\xi\cdot\mathbf{z})}.\tag{21}$$

Ψˆ <sup>+</sup> *<sup>μ</sup>* (**r**) are identified as the fields that are radiated forward into *z* > 0 by *ψ*ˆ*μ*(**x**). In (21), ˆ *ψ*˜*μ*(*ξ*) = ˆ *<sup>ψ</sup>*˜(*<sup>ξ</sup>* <sup>−</sup> *<sup>ξ</sup>***n**)*e*−*ikξ*·**xm** is the PW spectrum (5) of *<sup>ψ</sup>*ˆ*μ*(**x**), with <sup>ˆ</sup> *ψ*˜(*ξ*) being the spectrum of the "mother window" *ψ*ˆ(**x**). If *ψ*ˆ(**x**) is wide on a wavelength scale, then Ψˆ <sup>+</sup> *<sup>μ</sup>* (**r**) behave like collimated beams, emerging from the points **xm** over the *z* = 0 plane and directions ◦ *<sup>κ</sup>***<sup>n</sup>** = (*ξ***n**, *<sup>ζ</sup>***n**)=(sin *<sup>θ</sup>***<sup>n</sup>** cos *<sup>φ</sup>***n**, sin *<sup>θ</sup>***<sup>n</sup>** sin *<sup>φ</sup>***n**, cos *<sup>θ</sup>***n**) with *<sup>ζ</sup>***<sup>n</sup>** <sup>=</sup> <sup>1</sup> − |*ξ***n**|2.

#### *3.2. UWB Considerations*

In general, the applications in hand require UWB excitations. Following [17], we use the following frequency-scaling of the WFT frame set that renders the theory amenable for UWB field representations:

*(1) Frequency independent beam skeleton:* As implied from Equation (16) above, the beam lattice should be recalculated for each frequency. For efficient UWB representations, it is required to have the same beam lattice **xm**, *ξ***<sup>n</sup>** over the entire frequency band. In view of (16), this requirement implies (Equation (10) [21])

$$\nu(\omega) = \nu\_{\text{max}} \frac{\omega}{\omega\_{\text{max}}}, \qquad \omega \in [\omega\_{\text{min}}, \omega\_{\text{max}}]\_{\text{\textquotedblleft}\omega} \tag{22}$$

with *ν*max being the value of *ν* at *ω*max, so that *ν* < *ν*max for all *ω* < *ω*max. Typically, we use *ν*max = 1/3 (see discussion in (25) below).

*(2) Iso-diffracting propagators:* We use iso-diffracting (ID) Gaussian windows which are scaled with frequency in the form (Equation (27) [17])

$$
\hat{\mathfrak{gl}}\_{\rm ID}(\mathbf{x}) = e^{-k|\mathbf{x}|^{2}/2b}, \quad k > 0,\tag{23}
$$

where *b* > 0 is a frequency independent parameter. Inserting (23) into (21) and evaluating the integral one finds that the resulting propagators are ID-GBs, with *b* being the collimation distance. The ID designation of these Gaussian beams is due to the fact that their collimation distance *b* is frequency independent. This property implies that the beam propagation parameters are frequency independent even in inhomogeneous medium. Furthermore, when transformed into the TD, they give rise to ID-Pulsed beams (ID-PB) which are space time wave-packets that maintain their wave-packet structure even through propagation in inhomogeneous medium [49]. Explicit expressions for the corresponding phase–space beam propagators of (26) in free space are given in (Equations (28)–(29)) [17].

Typically *b* is chosen by the molder and depends on the application (see discussion in the numerical example below), but also should satisfy the condition *k*min*b* 1, which implies that the beams are highly collimated over the entire frequency band.

*(3) Snug frame:* The frame is optimal (or snug) when the frame elements are matched to the phase–space lattice (*x*¯, ¯ *ξ*) (in the sense that they should provide a balanced phase–space coverage). This requirement implies the relation *b* = *x*¯/ ¯ *ξ* [17]. Combining this condition with (16) one obtains (Equation (A2) [21]),

$$\left(\mathfrak{x},\mathfrak{f}\right) = \sqrt{\frac{2\pi\nu\_{\text{max}}}{k\_{\text{max}}}} \left(b^{1/2}, b^{-1/2}\right). \tag{24}$$

*(4) Simple expression for the dual frame function:* In view of (19) we have for *ν*max = 1/3 (Equation (A3) [21]),

$$\phi\_{\rm ID}(\mathbf{x}) \simeq \frac{\nu\_{\rm max}^2}{\pi b k\_{\rm max}^2} k^3 \,\hat{\psi}\_{\rm ID}(\mathbf{x}), \quad \omega \,\,\,\,\omega\_{\rm max}.\tag{25}$$

Over this regime *ϕ*ˆ is spatially and spectrally localized, and leads to a stable and localized expansion coefficients [17].

The properties above yield an efficient multi-frequency representation where the phase–space lattice and propagation parameters should be calculated only once for all frequencies in the band. These advantages also allow a simple transformation of the beam representation to the TD [19,21].

#### *3.3. The Beam Frame Theorem*

Following (21), we define the set of forward and backward propagators {Ψ<sup>ˆ</sup> <sup>±</sup> *<sup>μ</sup>* (**r**)} (compare Equation (21))

$$\Psi\_{\mu}^{\pm}(\mathbf{r}) = \left(\frac{k}{2\pi}\right)^{2} \int\_{|\xi|<1} d^{2}\xi \,\hat{\Psi}\_{\mu}(\xi) e^{ik(\xi \cdot \mathbf{x} \pm \zeta z)}, \quad |\xi\_{\mathbf{n}}| \le \xi\_{0} < 1. \tag{26}$$

where the parameter *ξ*<sup>0</sup> is typically chosen close to 1. Note that this subset includes only "propagating beams" whose spectrum, which is localized around *ξ***n**, is located in the propagating spectrum range |*ξ*| < 1. We denote this subset by the index *μP*. Inserting the ID Gaussian windows of (23) into (26) and evaluating the integrals asymptotically one readily identifies Ψˆ <sup>±</sup> *<sup>μ</sup>* as forward and backward ID-GB that propagate from *z* = ∓∞ to ±∞ in the directions ◦ *κ*± **<sup>n</sup>** = (*ξ***n**, ±*ζ***n**)=(sin *θ***<sup>n</sup>** cos *φ***n**, sin *θ***<sup>n</sup>** sin *φ***n**, cos *θ***n**) (see Figure 3), while for *z* = 0 they converge to the PS lattice of Section 3.1 as illustrated in Figure 3.

**Figure 3.** The forward/backward propagating beam frames Ψ± *<sup>μ</sup>* . (**a**) The forward and (**b**) the backward beam frames Ψ± *<sup>μ</sup>* . The BFs are illustrated by the hatched arrows. The ellipses correspond to the pulsed-beam-frames that are used in the TD formulations and are not considered here (see [21]).

As has been established by **the beam frame theorem** in [20], the beam-sets Ψˆ <sup>±</sup> *<sup>μ</sup>* (**r**) *μP* constitute frames (hence referred to as "beam frames" (BF)) at any *z* = *const*. plane over the Hilbert space H*<sup>P</sup>* of functions whose spectrum is bounded in the propagation domain <sup>|</sup>*ξ*<sup>|</sup> <sup>&</sup>lt; *<sup>ξ</sup>*0, with the set {Φ<sup>ˆ</sup> <sup>±</sup> *<sup>μ</sup>* (**r**)} being the dual frames. The propagators <sup>Φ</sup><sup>ˆ</sup> <sup>±</sup> *<sup>μ</sup>* have the same form as Ψˆ <sup>±</sup> *<sup>μ</sup>* in (26), except that <sup>ˆ</sup> *ψ*˜*<sup>μ</sup>* are now replaced by ˆ*ϕ*˜*μ*. Note that in view of (25), Φˆ <sup>±</sup> *μ* are proportional to Ψˆ <sup>±</sup> *μ* .

It follows from the beam frame theorem that any function over H*<sup>P</sup>* may be expanded by the BF. This observation is very important in the context of inverse scattering since it implies that both the scattered field and the medium are expanded on the same basis.

An important special case of the above is when the BF are used to expand forward or backward propagating wave-fields *u*ˆ±(**r**). In view of the theorem, *u*<sup>+</sup> may be expanded using Ψˆ <sup>+</sup> *<sup>μ</sup>* , and *u*<sup>−</sup> may be expanded using Ψˆ <sup>−</sup> *<sup>μ</sup>* , but the physically meaningful choice is to expand *u*<sup>±</sup> using Ψˆ <sup>±</sup> *<sup>μ</sup>* , respectively, viz (Equation (32) [20])

$$\mathcal{H}^{\pm}(\mathbf{r}) = \sum\_{\mu \in \mu \underline{\nu}} A\_{\mu}^{\pm} \Psi\_{\mu}^{\pm}(\mathbf{r}),\tag{27}$$

where the summation includes only "*μ<sup>P</sup>* propagating" frame-beams with no evanescent spectrum. As has been established by **the expansion coefficient invariance theorem** in [20], *A*ˆ<sup>±</sup> *<sup>μ</sup>* may be calculated by projecting *u*ˆ±(**r**) on the dual frame Φˆ <sup>±</sup> *<sup>μ</sup>* (**r**) over any *z* = *z* plane, giving the same result, i.e., (Equation (33) [20])

$$\mathcal{A}\_{\mu}^{\pm} = \left\langle \hat{\mathfrak{u}}^{\pm}(\mathbf{r}), \Phi\_{\mu}^{\pm}(\mathbf{r}) \right\rangle \Big|\_{z'} = \left\langle \hat{\mathfrak{u}}^{\pm}(\mathbf{x}), \Phi\_{\mu}^{\pm}(\mathbf{x}) \right\rangle \Big|\_{z=0} \tag{28}$$

where the last expression describes the canonical WFT coefficients of (18) evaluated over the *z* = 0 plane.

Finally we note that in [21], the BF theorem has been extended to the TD using ID-PB propagators.

#### **4. UWB Beam-Based Diffraction Tomography: Multi-Frequency Formulation**

The beam frame concept provides a framework to formulate the beam-based inverse scattering algorithm. Using the BFs, we may use the same set of beam basis functions to expand both the scattering data and the medium (actually, the sources that are induced due to the medium heterogeneities). As illustrated in Figure 4, the inverse problem is thereby described by the local interaction between the beam amplitudes and the unknown object. As noted in the introduction, optimal localization is obtained in the time-domain formulation, using localized space–time wave-packets. This, however, requires somewhat sophisticated processing tools [21]. In the present section we present only the multi-frequency formulation that utilizes conventional FD data-processing tools followed by integration over all the frequencies. The readers are referred to [36,37] for the TD interpretation.

**Figure 4.** The scattering mechanism within the propagating frame formulation. The incident field that propagates through the medium (see subplot (**a**)) gives rise to induced sources. At each *z* = *const*. plane, these sources are expanded by the forward/backward propagating BFs, giving rise to the forward/backward scattered fields depicted in subplot (**b**) in blue and red, respectively.

#### *4.1. The Inversion Algorithm*

Given the scattering data over the *zj* planes, the BF representation of the scattering fields into the *z* ≶ *zj* half spaces are given by (see (27))

$$\hat{\mu}\_{\rangle}^{s}(\mathbf{r},\omega) = \sum\_{\mu \in \mu p} A\_{\mu}^{\hat{j}}(\omega) \Psi\_{\mu}^{\mp}(\mathbf{r},\omega), \quad z \lessapprox' z\_{\rangle'} \tag{29}$$

where, as before, upper and lower signs correspond to *j* = 1 and *j* = 2, respectively. The expansion coefficients calculated via (28),

$$
\hat{A}\_{\mu}^{\hat{j}}(\omega) = \left< \hat{u}\_{\hat{j}}^{s}(\mathbf{x}, \omega), \Phi\_{\mu}^{\mp}(\mathbf{r}, \omega)|\_{z\_{\hat{j}}} \right>. \tag{30}
$$

Following the discussion after (7), these coefficients extract the local PW spectrum of the scattering data. Note that, as was done in the PW spectrum of Equation (5), the scattering WFT operation normalizes the scattering on the *zj* planes to their phase centers on the *z* = 0 plane. The coefficients in (30) are referred to as the beam-domain data.

The backpropagated fields *u*ˆ*<sup>b</sup> <sup>j</sup>*(**r**, *ω*) are obtained by extending (29) as is to *z* ≷ *zj* (see (II-9)). The "imaging fields" are then calculated by inserting (29) into Equation (11). In view of the local structure of the Ψˆ <sup>∓</sup> *<sup>μ</sup>* propagators, we obtain the explicit expression (II-11)

$$\mathcal{I}\_{\dot{\jmath}}(\mathbf{r}) \simeq \frac{2}{i\omega} e^{-ik\mathbf{\hat{x}}^{\ddot{\jmath}}\mathbf{\hat{r}}} \sum\_{\mu} A\_{\mu}^{\dot{\jmath}}(\omega) \cos^{2}\left(\frac{\gamma\_{\mathbf{n}}^{\ddagger}}{2}\right) \Psi\_{\mu}^{\ddagger}(\mathbf{r}, \omega), \tag{31}$$

where *γ*∓ **<sup>n</sup>** represents the angle between the illumination direction <sup>−</sup>◦ *κ<sup>i</sup>* and the beam direction ◦ *κj <sup>μ</sup>* (which actually depends only on **n**. Finally, the reconstructed object is calculated via (12) and (13). For full details, the reader should refer to Appendices II-A,B.

#### *4.2. Interpretation within the Born Approximation*

In order to gain insight into the structure of the beam-domain data, we insert the Born approximation of the scattered field in (7) into (30). The resulting relation between *O*(**r**) and the beam data is given by (I-21)

$$A^j\_{\mu}(\omega) \simeq \langle \mathcal{O}(\mathbf{r}), \,\, \hat{\Lambda}^j\_{\mu}(\mathbf{r}, \omega) \rangle\_{\mathbb{V}'}, \quad \hat{\Lambda}^j\_{\mu}(\mathbf{r}, \omega) = k/2i \cos \theta\_{\mathbf{n}} e^{-i \mathbf{k} \mathbf{\bar{x}'} \mathbf{r}} \Phi^{\mp}\_{\mu}(\mathbf{r}, \omega), \tag{32}$$

where the integration covers the entire scatterer domain. Thus, within the Born approximation, the data is described as projections of *O*(**r**) on the beam axis, using the projection kernels Λˆ *<sup>j</sup> <sup>μ</sup>*(**r**, *ω*). As shown in I-Section 5-B in [36], this projection extracts the local stratification of *O* along the beam axis at an angle *γ*∓ **<sup>n</sup>** defined in (31). This implies that the scattering amplitudes *A*ˆ*<sup>j</sup> <sup>μ</sup>* are obtained from Snell's reflections due the stratification components in *O*(**r**) along the *μ* beam axis, so that strong amplitudes are obtained only for those *μ* (locations and directions) that correspond to the stratification of *O*(**r**) along the *μ* beam axis. Note that (32) is the local generalization of (7), where the BF basis is used instead of the conventional Green's function that radiates in all directions.

Further localization along the beam axis is provided by using the TD formulation in (I-27)-(I-32). However, as noted earlier, the TD formulation is not presented here since our goal in this paper is to present the pragmatic and practical formulation in the FD where all the operations are based on Fourier-type integrals. The readers are referred to [36,37] for more details on the TD formulations.

To further explore the FD beam data representation, we consider the spectral representation of *A*ˆ*<sup>j</sup> <sup>μ</sup>*. Substituting (5) into (30) and changing the order of integrations, *<sup>A</sup>*ˆ*<sup>j</sup> <sup>μ</sup>* can be expressed as (I-22)

$$A^j\_{\mu} \simeq \left(\frac{k}{2\pi}\right)^2 \int d^2\xi \, \frac{ik}{2\mathcal{J}} \hat{\phi}^\*\_{\mu}(\xi) \mathcal{O}(\mathbf{K})\Big|\_{\mathbf{K}=k(\stackrel{\circ}{\mathbf{k}}\_j-\stackrel{\circ}{\mathbf{k}}^j)}.\tag{33}$$

The expression is the local alternative to the plane-wave-based DT identity in (8). It is recognized as ◦ *κ* ± **<sup>n</sup>** samples of the value of *O*¯(**K**) over the shifted Ewald sphere defined in (9). The spectral width of these samples is that of ˆ*ϕ*˜*μ*(*ξ*).

#### **5. Numerical Examples**

In this section, we demonstrate the beam-based DT algorithm through a numerical example. We begin with a step-by-step summary of the algorithm, including guidelines for choosing the parameters.

#### *5.1. A Step by Step Summary of the Algorithm*

#### **Phase I—the experimental setup**

1. We consider a realistic case where the object is illuminated by an array of *M* independent point transducers over the *z*<sup>1</sup> plane, as illustrated in Figure 5a. We illustrate here only the reflected data on the *z*<sup>1</sup> plane, since in many applications (e.g., geophysics) the transmission field cannot be measures, and in some cases this data is not needed (see discussion in Section 2.4). If, however, the transmitted data at *z*<sup>2</sup> is available, then the receiver array considerations are similar.


One may also calculate the time-harmonic PW spectrum of the scattered field via *<sup>p</sup>*-stacking with proper phase terms. The result is an *<sup>M</sup>* <sup>×</sup> *<sup>M</sup>* data matrix <sup>ˆ</sup> **U**˜ *s <sup>p</sup>q*(*ω*) describing the *p* spectral PW due to an excitation by the *q* incident PW. As noted earlier, before we do the stacking, the array dimensions should be zero- padded to avoid aliasing of the final image when the images are generated by spectral integrations.


**Figure 5.** Guidelines for choosing the experimental setup. (**a**) The size of the array should be wide enough to provide a plane-wave illumination at the target range *R*. The scan angle is limited in order to be sufficiently far for the end-point diffraction zone. (**b**) For local imaging, only the part of the object inside the DOI should satisfy the conditions above.

#### **Phase II—Constructing the phase space lattice**

The next step is to set up the phase–space lattice and choose the expansion parameters


**Phase III—Calculating the beam data**


**Phase IV—Local reconstruction via beam backpropagation**


**Figure 6.** An overview of the local inversion algorithm. The beam expansion of the scattered field is plotted as gray arrows. Only those covering the array are plotted. The scattering data is then transformed to beam amplitudes by stacking the receivers data via (30), as schematized by the black ellipses. Only beams with high amplitude are considered and backpropagated via (31). The image in the DoI (black rectangle) is obtained by aggregating the contribution of beams that pass inside the DoI (red beams).

#### *5.2. Example A: A Smooth and Quasi Stratified Medium. UWB Reflection Mode Data*

The medium is plotted in Figure 7a in a 2*D* coordinate frame **r** = (*x*, *z*), with the DoI being the 20 × 20 black rectangle. For simplicity we normalize all space-units such that the background wave-speed is *v*<sup>0</sup> = 1. Note that the contrast is relatively large with values of *O*max = ±0.44. Note that this example is one of those treated in [37] (see Figure 6, but here the processing is done in the multi-frequency domain as outlined above.

The medium is dominated by stratification along the *z* direction, hence its *K*-space distribution is localized near the *Kz* axis (see discussion below). Referring to the discussion in Section 2.4, it can be recovered using UWB reflection data on the *z*<sup>1</sup> plane. We therefore use illumination by a *z*-propagating time-harmonic PW with frequencies in the band Ω = 0.1, 1 . The exact data is generated using method of moments (MoM) simulations. We record only the reflected data over an array of receivers located at *z* = −150 with |*x*| < 250 with inter-element spacing *d* = 1.15 *π* (larger than the Nyquist distance).

We set *b* = 50, such that the beams are collimated inside the DoI, while maintaining *k*min*b* 1 as required for collimation after (23). Using also *k*max = 1 and *ν*max = 1/3 we obtain from (24) (*x*¯, ¯ *ξ*)=(9.71, 0.194). The resulting BF propagators Ψˆ <sup>±</sup> *<sup>μ</sup>* (**r**), Φˆ <sup>±</sup> *<sup>μ</sup>* (**r**) are calculated via (Equations (C1)–(C5) [21]).

Next we calculate the beam-domain data *A*ˆ*<sup>j</sup> <sup>μ</sup>* via (30). The reconstructed object inside the DoI is found using the reflected data imaging field ˆ*I*1(*ω*) of (31) where we consider only backpropagated beams whose *μ* axis passes inside the DoI, and then integrating over all frequencies as in (12). The reconstructed medium is illustrated in Figure 7b. As can be seen, the reconstructed object matches well with the object inside the DoI. To better quantify the image results, in Figure 8 we plot cross-sectional cuts of the object at *x* = 0 and at *x* ± 6.

**Figure 7.** Example A: (**a**) The original (**a**) and the reconstructed (**b**) object functions. The DoI is illustrated by the black rectangle.

**Figure 8.** Example A: Comparison of the results along cross sectional cuts parallel to the *z*-axis.

The sources of error are readily seen in the *K*−space distribution of the original and reconstructed media in Figure 9. Note that the imaging algorithm has recovered most of the object's *K*−space, except for the region around |**K**| ≈ 0. As discussed in Section 2.4, this missing data is due to the low frequency cutoff *kmin* = 0.1 in the data. The main drawback of the "UWB reflection mode inversion" schemes is the missing transversal spectrum components and the |**K**| → 0 components, (which are small in this example). In general, one may try to recover this data by using transition mode data (*z*2) but in many applications this data is not available. Alternatively, one may use several illumination directions which are synthesized from the array data via the method outlined in Phase I of Section 5. These additional illuminations are also used to reduce the reconstruction noise, as explored in II-Section 6 in [37]. The readers are referred to other examples in [36,37,42–44,50].

**Figure 9.** Example A: Comparison of the *K*− space distributions of the original (**a**) and reconstructed (**b**) objects.

#### *5.3. Example B: An Object with Sharp Boundaries. Reflection and Transmission Data*

The object shown in Figure 10a has sharp boundaries, strong transversal *K* components, and a non-zero average (i.e., *<sup>O</sup>*¯(|**K**<sup>|</sup> <sup>=</sup> <sup>0</sup>) <sup>=</sup> 0). As before, we consider a 2*<sup>D</sup>* problem with **r** = (*x*, *z*) and normalize the space-units such that the background wave-speed is *v*<sup>0</sup> = 1.

The source array is located on the *z*<sup>1</sup> = −150 plane over |*x*| < 250, with interelement spacing *d* = 1.15*π*. Using this array, we may synthesize PW illumination over a spectral range of ±60◦ (see discussion in Section 5, Phase I(1–4)). The frequency band is Ω = 0.1, 1 . We consider both the reflection and transmission data over similar receiver arrays at *z*<sup>1</sup> = −150 and *z*<sup>2</sup> = 150. The exact scattered data is calculated via the MoM.

For the beam processing we use *b* = 20, such that the beams are collimated inside the DoI, while maintaining *k*min*b* 1 as required for collimation after (23). Using also *k*max = 1 and *ν*max = 1/3 we obtain from (24) (*x*¯, ¯ *ξ*)=(6.47, 0.32). The resulting BF propagators Ψˆ <sup>±</sup> *<sup>μ</sup>* (**r**), Φˆ <sup>±</sup> *<sup>μ</sup>* (**r**) are calculated via (Equations (C1)–(C5) [21]).

Figure 10b depicts the reconstructed objects in the front (left) using a single PW illumination at *θ<sup>i</sup>* = 0 and reflection data at the *z*<sup>1</sup> plane. As expected, the reflection data provides good longitudinal resolution but poor transversal resolution (see Figure 2b). Note also that the value of the reconstructed object function is approximately one half of the true value due to the missing data at |**K**| = 0. As expected, the reconstruction of the object outside the DoI is poor.

In Figure 10c we improved the resolution by using several illumination directions (as one may expect by considering Figure 2a,b), yet the reconstruction still suffers from poor transversal resolution and low value of the reconstructed object. These problems are mitigated in Figure 10d where we used both the reflection and transmission data as in (13). Further improvement can be made via iterative schemes [36,37,42–44,50].

Finally, in Figure 11 we demonstrate local imaging within different DoIs. Figure 11a depicts the reconstruction of a cylinder at the front (left-top) using both reflection and transmission data due to illumination at *θ<sup>i</sup>* = 0. The reconstruction is a bit poorer than in Figure 10d where we used several illumination directions. As expected, the reconstruction outside this DoI is poor.

**Figure 10.** Numerical example B. (**a**) The object function. (**b**) Reconstruction using single illumination and reflection data. (**c**) Reconstruction using several illuminations and reflection data. (**d**) Reconstruction using several illuminations and reflection and transmission data. The DoI is illustrated in a black rectangle.

Figure 11b depicts the reconstruction of a cylinder at the back (right-top). Since this cylinder is poorly illuminated by normal incidence, we use here reflection and transmission data from several illumination directions *<sup>θ</sup><sup>i</sup>* <sup>=</sup> <sup>∓</sup>30◦, <sup>∓</sup>40◦. The reconstruction inside the DoI is much better than in Figure 10. As before, the reconstruction outside this DoI is poor.

**Figure 11.** Numerical example B: Local reconstruction in different DoI's (black rectangles). (**a**) Reconstruction of a cylinder at the front using both reflection and transmission data due to illumination at *θ<sup>i</sup>* = 0. (**b**) Reconstruction of a cylinder in the back using both reflection and transmission data from several illumination directions.

#### **6. Discussion and Conclusions**

In this paper, we reviewed the local diffraction–tomography inversion scheme introduced originally in [36,37]. The method is based on a local transformation of the scattering data and local reconstruction using beam backpropagation. It is structured on the concept of beam-frames (BFs). The BFs are a phase–space set of beam-waves that constitute local basis functions (frames) over the propagation domain. As such, they provide an alternative local basis for the global PW or Green's function radiation integrals. We use the BFs to formulate a local inversion algorithm as an alternative to the conventional approaches. In this and other publications, we demonstrated and explored the advantages of the local algorithm over the conventional PW and Green's function DT algorithms:


In previous publications [36,37] we have emphasized TD data processing, where the beam waves are localized space–time wave-packets. This requires somewhat sophisticated mathematics to construct the wave-packets and use them as signal processing tools. The main advantage of the TD approach is the data localization and interpretation. In the present paper, on the other hand, we utilized FD processing followed by an integration over the relevant frequency band. The motivation has been to provide the reader with more straightforward Fourier-based data-processing tools. We also provide the processing tools and closed-form expressions for the local imaging formula, as well as step-by-step guidelines for choosing the various scheme parameters. The method provides an efficient UWB formulation where one has to calculate the beam lattice and propagators only once and then use them for all the frequencies.

**Funding:** This work is partially supported by the Israeli Science Foundation (ISF) under Grant Grant 1111/19.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data associated with this research are available and can be obtained by contacting the author.

**Acknowledgments:** The author would like to express his sincere appreciation to E. Heyman for helpful discussions pertaining to the preparation of this manuscript.

**Conflicts of Interest:** The author declare no conflict of interest.

#### **References**

