A Multi-Frequency Tomographic Inverse Scattering Using Beam Basis Functions

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 multi-frequency 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.

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,2,3,4,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 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,14,15,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,23,24,25,26,27,28,29,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,32,33,34,35,36,37] and [38,39,40,41,42,43,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).

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:

a. 

Data localization: The phase–space processing of the scattered data extract the local direction of arrival. The phase space representation of the data is therefore localized along well defined trajectories corresponding to the local direction and time of arrival from the relevant regions in the target domain.

b. 

Under the Born approximation, the beam-wave scattering mechanism by the medium inhomogeneities is described by local Snell’s reflections from the local stratification, which is related to the local Radon transform (LRT) of the medium inhomogeneities (Section 5 in [36]).

c. 

As follows from the discussion in items a and b, the phase space data is directly related to the LRT of the medium inhomogeneities about a given region.

d. 

The beam-based imaging enables local imaging within a given domain of interest (DoI) by considering only the data that correspond to beams that pass through or near the DoI. This not only reduces the problem complexity, but also reduces the noise level, since data and noise arriving from other regions are a priori filtered out.

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 $\mathring{\mathbf{\kappa }}{}^i$ 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 $\mathring{\mathbf{\kappa }}{}^i$ (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_1<0$ and at $z=z_2>0$. We assume a $3D$ coordinate frame $\mathbf{r}=(\mathbf{x},z)$ where the z-coordinate is normal to the measurement planes, and $\mathbf{x}=(x_1;x_2)$ are the transversal coordinates. The data is collected over a wide frequency band $\Omega \in [{\omega }_{\min },{\omega }_{\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}\left(\omega \right)=\int dt u(t)e^{i\omega t},$$

(1)

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

The unknown object is embedded in a uniform background wavespeed $v_0$ and assumed to be lossless and non-dispersive. It is described by the unknown wavespeed $v\left(\mathbf{r}\right)$, and we define the “object function”

$$O\left(\mathbf{r}\right)=(v_0/v(\mathbf{r}){)}^2-1$$

(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

$$\hat{u}{}^i(\mathbf{r},\omega )=\hat{F}(\omega )e^{ik{\mathring{\mathbf{\kappa }}}^i·\mathbf{r}}, u^i(\mathbf{r},t)=F(t-v_0^{-1}{\mathring{\mathbf{\kappa }}}^i·\mathbf{r}),$$

(3)

where $\hat{F}(\omega )$ is the source spectrum and $k=\omega /v_0$ is the wavenumber. The incident wave propagates in the direction

$${\mathring{\mathbf{\kappa }}}^i=({\mathbf{\xi }}^i,{\zeta }^i)=\sin {\theta }^i \cos {\varphi }^i{\mathring{x}}_1+\sin {\theta }^i \sin {\varphi }^i{\mathring{x}}_2+\cos {\theta }^i\mathring{z},$$

(4)

with $({\theta }^i,{\varphi }^i)$ being the polar angles with respect to the z axis, and over-circles denote unit vectors.

The scattered fields measured over the $z_j$ planes, $j=1,2$ are denoted as ${\hat{u}}_j^s(\mathbf{x},\omega )$ (see Figure 1a). The PW spectral representation of ${\hat{u}}_j^s\left(\mathbf{r}\right)$ is defined via

$$\hat{\tilde{u}}{}_j^s(\mathbf{\xi })=e^{\pm ik\zeta z_j}\int d^{2}x\hat{u}{}_j^s(\mathbf{x})e^{-ik\mathbf{\xi }·\mathbf{x}}, \zeta =\sqrt{1-\mathbf{\xi }·\mathbf{\xi }}, \mathrm{Im}\zeta \ge 0.$$

(5)

where lower and upper signs correspond to $j=1$ and 2, respectfully. We added the $e^{\pm ik\zeta z_j}$ phase term in (5) in order to normalize the spectral PWs to the $z=0$ plane instead of $z=z_j$ planes.

Note that we use here the frequency normalized spectral coordinates $\mathbf{\xi }={\mathbf{K}}_x/k$ which are related to the PW direction via $\mathbf{\xi }=({\xi }_1,{\xi }_2)=\sin \theta (\cos \varphi ,\sin \varphi )$, where $(\theta ,\varphi )$ are the conventional spherical angles with respect to the z-axis so that the scattered PWs propagate in the unit vector direction

$${\mathring{\mathbf{\kappa }}}_j=(\mathbf{\xi },\mp \zeta )=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta ), j=1,2.$$

(6)

The spectral ranges $\left|\mathbf{\xi }\right|<1$ and $\left|\mathbf{\xi }\right|>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}^3{r}^{\prime }{\hat{u}}^i({\mathbf{r}}^{\prime })O({\mathbf{r}}^{\prime })\hat{G}(\mathbf{r},{\mathbf{r}}^{\prime }),$$

(7)

where $\hat{G}=\frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}$ is the 3D Green’s function in the uniform background. This approximation is valid if $O\left(\mathbf{r}\right)\ll 1$ and in addition $kL{O}_{\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 $\hat{G}$, we obtain (I-7),

$$\hat{\tilde{u}}{}_j^s(\mathbf{\xi })\simeq \frac{k}{-2i\zeta }\overline{O}(\mathbf{K}){\mid }_{\mathbf{K}=k({\mathring{\mathbf{\kappa }}}_j-{\mathring{\mathbf{\kappa }}}^i)} |\mathbf{\xi }|<1,$$

(8)

where ${\mathring{\mathbf{\kappa }}}_j$ and ${\mathring{\mathbf{\kappa }}}^i$ are given by (6) and (4), and

$$\overline{O}(\mathbf{K})=\int d^{3}rO(\mathbf{r})e^{-i\mathbf{K}·\mathbf{r}}, \mathbf{K}=(K_1,K_2,K_z),$$

(9)

is the K-space distribution of $O\left(\mathbf{r}\right)$. Equation (8) is referred to as the DT identity. It relates the scattering data in the ${\mathring{\mathbf{\kappa }}}_j$ directions to values of $\overline{O}\left(\mathbf{K}\right)$ at the points $\mathbf{K}=k({\mathring{\mathbf{\kappa }}}_j-{\mathring{\mathbf{\kappa }}}^i)$. As illustrated in Figure 1b, these points define a K-space sphere with radius k that is centered at $\mathbf{K}=-k{\mathring{\mathbf{\kappa }}}^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 $z_j$ measurement plane with $j=1,2$, respectively.

The DT identity above applies only to the propagation spectrum $\left|\mathbf{\xi }\right|<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 ${\mathring{\mathbf{\kappa }}}^i$ 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 $\overline{O}\left(\mathbf{K}\right)$. Note that for lossless (real) objects, $\overline{O}\left(\mathbf{K}\right)={\overline{O}}^{*}(-\mathbf{K})$, so that only half of the K-space needs to be recovered.

As one observes from Figure 2a, the transmitted data on $z_2$ recovers the K-space distribution of the object in a sphere of radius $\sqrt{2}k$ 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 $\overline{O}\left(\mathbf{K}\right)$. Note that in the limit of $k\to \infty$, the transmitted data hemispheres in Figure 2a reduce to planar surfaces normal to ${\mathring{\mathbf{\kappa }}}^i$ that pass through the origin, thus providing the K-space representation of conventional X-ray tomography [46].

One option to reconstruct $O\left(\mathbf{r}\right)$ is to recover $\overline{O}\left(\mathbf{K}\right)$ 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 $\overline{O}\left(\mathbf{K}\right)$ 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 $\Omega \in [{\omega }_{\min },{\omega }_{\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 $\pi /4$ cone with an axis along the negative $K_z$ 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_1$ plane mainly recovers the object variations along the z axis, while the transmitted data on the $z_2$ 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_1$ plane, but may not be sufficient for objects with a large transversal variations. Another limitation is the missing data for $|K_z|<2k_{\min }$ (see Figure 2b), while the missing data for $|K_z|>2k_{\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 $\left|\mathbf{K}\right|$. Note also that the maximal axial resolution for the case of normal incidence is ${\delta }_z=\pi /k_{\max }$.

The object can be reconstructed using an inverse transform of $\overline{O}\left(\mathbf{K}\right)$. 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 ${\hat{u}}^s(\mathbf{x},\omega )$ over the $z_j$ planes, the backpropagated fields into the $z>z_1$ and $z<z_2$ regions are given by (see (5))

$$\hat{u}{}_j^b(\mathbf{r})=(\frac{k}{2\pi }{)}^2{\int }_{\left|\xi \right|<1}{d}^2\xi \hat{\tilde{u}}{}_j^s(\mathbf{\xi })e^{ik(\mathbf{\xi }·\mathbf{x}\mp \zeta z)},$$

(10)

where we restrict the integration to the visible spectrum $\left|\mathbf{\xi }\right|<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}{}_j(\mathbf{r},\omega )=v_0^{-1}{k}^{-2}{\mathring{\mathbf{\kappa }}}^i·\nabla [e^{-ik{\mathring{\mathbf{\kappa }}}^i·\mathbf{r}}\hat{u}{}_j^b(\mathbf{r},\omega )].$$

(11)

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

$$\breve{O}(\mathbf{r})=2\mathrm{Re}\frac{1}{\pi }{\int }_0^{\infty }d\omega \hat{I}{}_j(\mathbf{r},\omega ).$$

(12)

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

$$\breve{O}(\mathbf{r})={\breve{O}}_1(\mathbf{r})+{\breve{O}}_2(\mathbf{r}).$$

(13)

The features of $O\left(\mathbf{r}\right)$ that are described by ${\breve{O}}_j$ have been discussed above in connection with Figure 2b. As noted there, in many situations it is sufficient to recover only ${\breve{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 $\mathbf{r}=(\mathbf{x},z)$, $\mathbf{x}=(x_1,x_2)$, due to a time harmonic field ${\hat{u}}_0(\mathbf{x},\omega )$ defined over the plane $z=0$.

The WFT frame set $\left\{{\hat{\psi }}_{\mathit{\mu }}(\mathbf{x},\omega )\right\}$ is defined by (Equation (22)] [17])

$${\hat{\psi }}_{\mathit{\mu }}(\mathbf{x},\omega )=\hat{\psi }(\mathbf{x}-\mathbf{m}\overline{x})e^{ik\mathbf{n}\overline{\xi }·(\mathbf{x}-\mathbf{m}\overline{x})},$$

(14)

with $\hat{\psi }$ being a localized window function (typically a Gaussian, see more details below) and $\mathit{\mu }=\left(\mathbf{m},\mathbf{n}\right)$ being a 4-index. The frame elements are localized about the spatial $\left(\mathbf{x}\right)$ and spectral $\left(\mathbf{\xi }\right)$ phase space lattice

$$({\mathbf{x}}_{\mathbf{m}}, {\mathbf{\xi }}_{\mathbf{n}})=(\mathbf{m}\overline{x}, \mathbf{n}\overline{\xi })=(m_1\overline{x},m_2\overline{x}; n_1\overline{\xi },n_2\overline{\xi }),$$

(15)

where $(\overline{x},\overline{\xi })$ defines the lattice unit cell. As will be shown, the points $({\mathbf{x}}_{\mathbf{m}},{\mathbf{\xi }}_{\mathbf{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\pi$, implying that

$$k\overline{x}\overline{\xi }=2\pi \nu ,$$

(16)

with $\nu <1$ being the overcompleteness parameter and the limit $\nu =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 ${\hat{u}}_0\left(\mathbf{x}\right)$ in the form

$${\hat{u}}_0(\mathbf{x})=\sum _{\mathit{\mu }}{\hat{a}}_{\mathit{\mu }}{\hat{\psi }}_{\mathit{\mu }}(\mathbf{x}).$$

(17)

In view of the overcompleteness, the coefficients set $\left\{{\hat{a}}_{\mathit{\mu }}\right\}$ is not unique. A particularly convenient set with a minimum ${\ell }_2$ norm is obtained by using the dual frame $\left\{{\hat{\phi }}_{\mathit{\mu }}(\mathbf{x})\right\}$ which has the same structure as $\left\{{\hat{\psi }}_{\mathit{\mu }}\right\}$ in (15) except that the mother window $\hat{\psi }(\mathbf{x})$ is replaced by the dual mother window $\hat{\phi }\left(\mathbf{x}\right)$. The resulting canonical coefficient set is given by (Equation (23) [17]).

$${\hat{a}}_{\mathit{\mu }}=\langle {\hat{u}}_0(\mathbf{x}),{\hat{\phi }}_{\mathit{\mu }}(\mathbf{x})\rangle$$

(18)

where $〈f,g〉=\int f{g}^{*}$ is the conventional ${\mathbb{L}}_2$ inner product in the transverse coordinate x. The canonical coefficients ${\hat{a}}_{\mathit{\mu }}$ in (18) are readily identified as the local spectrum of ${\hat{u}}_0\left(\mathbf{x}\right)$ windowed with respect to ${\hat{\phi }}_{\mathit{\mu }}$ about the phase–space points $\left({\mathbf{x}}_{\mathbf{m}},{\mathbf{\xi }}_{\mathbf{n}}\right)$.

Generally, $\hat{\phi }$ should be calculated numerically, for a given $\hat{\psi }$ and lattice $\left(\overline{x},\overline{\mathbf{\xi }}\right)$. However, if the lattice is sufficiently overcomplete, ($\nu \lesssim 1/3$) $\hat{\phi }\propto \hat{\psi }$ can be approximated by (Equation (11) [17])

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

(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): $\left(i\right)$—as follows from (19), in this parameter regime $\hat{\phi }$ is localized both spatially and spectrally, so that the expansion (17) comprises local and stable coefficients. $\left(ii\right)$—$\hat{\phi }$ is given analytically via (19) and does not have to be to calculated numerically. Reason $\left(ii\right)$ is critical for UWB problems where $\hat{\phi }$ needs to be found for each $\omega$.

The radiated field in $z>0$ is obtained now by replacing ${\hat{\psi }}_{\mathit{\mu }}\left(\mathbf{x}\right)$ in Equation (17) by beam propagators (Equation (24) [17])

$$\hat{u}(\mathbf{r})=\sum _{\mathit{\mu }}{\hat{a}}_{\mathit{\mu }}{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{+}(\mathbf{r}),$$

(20)

$${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{+}(\mathbf{r})=(\frac{k}{2\pi }{)}^2\int d^2\xi \hat{\tilde{\psi }}{}_{\mathit{\mu }}(\mathbf{\xi })e^{ik(\mathbf{\xi }·\mathbf{x}+\zeta z)}.$$

(21)

${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{+}(\mathbf{r})$ are identified as the fields that are radiated forward into $z>0$ by ${\hat{\psi }}_{\mathit{\mu }}\left(\mathbf{x}\right)$. In (21), $\hat{\tilde{\psi }}{}_{\mathit{\mu }}(\mathbf{\xi })=\hat{\tilde{\psi }}(\mathbf{\xi }-{\mathbf{\xi }}_{\mathbf{n}})e^{-ik\mathbf{\xi }·{\mathbf{x}}_{\mathbf{m}}}$ is the PW spectrum (5) of ${\hat{\psi }}_{\mathit{\mu }}\left(\mathbf{x}\right)$, with $\hat{\tilde{\psi }}(\mathbf{\xi })$ being the spectrum of the “mother window” $\hat{\psi }\left(\mathbf{x}\right)$. If $\hat{\psi }\left(\mathbf{x}\right)$ is wide on a wavelength scale, then ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{+}(\mathbf{r})$ behave like collimated beams, emerging from the points ${\mathbf{x}}_{\mathbf{m}}$ over the $z=0$ plane and directions ${\mathring{\mathbf{\kappa }}}_{\mathbf{n}}=({\mathbf{\xi }}_{\mathbf{n}},{\zeta }_{\mathbf{n}})=(\sin {\theta }_{\mathbf{n}}\cos {\varphi }_{\mathbf{n}},\sin {\theta }_{\mathbf{n}}\sin {\varphi }_{\mathbf{n}},\cos {\mathrm{\theta }}_n)$ with ${\zeta }_{\mathbf{n}}=\sqrt{1-|{\mathbf{\xi }}_{\mathbf{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 $\left({\mathbf{x}}_{\mathbf{m}},{\mathbf{\xi }}_{\mathbf{n}}\right)$ over the entire frequency band. In view of (16), this requirement implies (Equation (10) [21])

$$\nu (\omega )={\nu }_{\max }\frac{\omega }{{\omega }_{\max }}, \omega \in [{\omega }_{\min },{\omega }_{\max }],$$

(22)

with ${\nu }_{\max }$ being the value of $\nu$ at ${\omega }_{\max }$, so that $\nu <{\nu }_{\max }$ for all $\omega <{\omega }_{\max }$. Typically, we use ${\nu }_{\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{\psi }}_{\mathrm{ID}}(x)=e^{-k|\mathbf{x}{|}^2/2b}, k<0,$$

(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\gg 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 $(\overline{x},\overline{\xi })$ (in the sense that they should provide a balanced phase–space coverage). This requirement implies the relation $b=\overline{x}/\overline{\xi }$ [17]. Combining this condition with (16) one obtains (Equation (A2) [21]),

$$(\overline{x},\overline{\xi })=\sqrt{\frac{2\pi {\nu }_{\max }}{k_{\max }}}(b^{1/2}, {-b}^{1/2}).$$

(24)

(4) Simple expression for the dual frame function: In view of (19) we have for ${\nu }_{\max }=1/3$ (Equation (A3)[21]),

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

(25)

Over this regime $\hat{\phi }$ 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 $\{{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }(\mathbf{r})\}$ (compare Equation (21))

$${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }(\mathbf{r})=(\frac{k}{2\pi }{)}^2{\int }_{\left|\mathbf{\xi }\right|<1}{d}^2\xi \hat{\tilde{\psi }}{}_{\mathit{\mu }}(\mathbf{\xi })e^{ik(\mathbf{\xi }·\mathbf{x}\pm \zeta z)}, |{\mathbf{\xi }}_{\mathbf{n}}|\le {\xi }_0<1.$$

(26)

where the parameter ${\xi }_0$ is typically chosen close to 1. Note that this subset includes only “propagating beams” whose spectrum, which is localized around ${\mathbf{\xi }}_{\mathbf{n}}$, is located in the propagating spectrum range $\left|\mathbf{\xi }\right|<1$. We denote this subset by the index ${\mathit{\mu }}_P$. Inserting the ID Gaussian windows of (23) into (26) and evaluating the integrals asymptotically one readily identifies ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }$ as forward and backward ID-GB that propagate from $z=\mp \infty$ to $\pm \infty$ in the directions $\mathring{\mathbf{\kappa }}{}_{\mathbf{n}}^{\pm }=({\mathbf{\xi }}_{\mathbf{n}},\pm {\zeta }_{\mathbf{n}})=(\sin {\theta }_{\mathbf{n}}\cos {\varphi }_{\mathbf{n}},\sin {\theta }_{\mathbf{n}}\sin {\varphi }_{\mathbf{n}},\cos {\theta }_{\mathbf{n}})$ (see Figure 3), while for $z=0$ they converge to the PS lattice of Section 3.1 as illustrated in Figure 3.

As has been established by the beam frame theorem in [20], the beam-sets ${\left\{{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right)\right\}}_{{\mathit{\mu }}_P}$ constitute frames (hence referred to as “beam frames” (BF)) at any $z=const.$ plane over the Hilbert space ${\mathbb{H}}_P$ of functions whose spectrum is bounded in the propagation domain $\left|\mathbf{\xi }\right|<{\xi }_0$, with the set $\left\{{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right)\right\}$ being the dual frames. The propagators ${\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }$ have the same form as ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }$ in (26), except that ${\hat{\tilde{\psi }}}_{\mathit{\mu }}$ are now replaced by ${\hat{\tilde{\phi }}}_{\mathit{\mu }}$. Note that in view of (25), ${\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }$ are proportional to ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }$.

It follows from the beam frame theorem that any function over ${\mathbb{H}}_P$ 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 ${\hat{u}}^{\pm }\left(\mathbf{r}\right)$. In view of the theorem, $u^{+}$ may be expanded using ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{+}$, and $u^{-}$ may be expanded using ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{-}$, but the physically meaningful choice is to expand $u^{\pm }$ using ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }$, respectively, viz (Equation (32) [20])

$${\hat{u}}^{\pm }(\mathbf{r})=\sum _{\mathit{\mu }\in {\mathit{\mu }}_P}{\hat{A}}_{\mathit{\mu }}^{\pm }{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }(\mathbf{r}),$$

(27)

where the summation includes only “${\mathit{\mu }}_P$ propagating” frame-beams with no evanescent spectrum. As has been established by the expansion coefficient invariance theorem in [20], ${\hat{A}}_{\mathit{\mu }}^{\pm }$ may be calculated by projecting ${\hat{u}}^{\pm }\left(r\right)$ on the dual frame ${\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }(\mathbf{r})$ over any $z=z^{\prime }$ plane, giving the same result, i.e., (Equation (33) [20])

$${\hat{A}}_{\mathit{\mu }}^{\pm }=\langle {\hat{u}}^{\pm }(\mathbf{r}),{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }(\mathbf{r})\rangle {\mid }_{z\prime }=\langle {\hat{u}}^{\pm }(\mathbf{x}),{\hat{\phi }}_{\mathit{\mu }}^{\pm }(\mathbf{x})\rangle {\mid }_{z=0}$$

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

4.1. The Inversion Algorithm

Given the scattering data over the $z_j$ planes, the BF representation of the scattering fields into the $z\lessgtr z_j$ half spaces are given by (see (27))

$${\hat{u}}_j^s(\mathbf{r},\omega )=\sum _{\mathit{\mu }\in {\mathit{\mu }}_P}{\hat{A}}_{\mathit{\mu }}^j(\omega ){\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\mp }(\mathbf{r},\omega ), z\lessgtr z_j,$$

(29)

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

$${\hat{A}}_{\mathit{\mu }}^j(\omega )=\langle {\hat{u}}_j^s(\mathbf{x},\omega ),{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\mp }(\mathbf{r},\omega ){\mid }_{z_j}\rangle .$$

(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 $z_j$ 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 ${\hat{u}}_j^b(\mathbf{r},\omega )$ are obtained by extending (29) as is to $z\gtrless z_j$ (see (II-9)). The “imaging fields” are then calculated by inserting (29) into Equation (11). In view of the local structure of the ${\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\mp }$ propagators, we obtain the explicit expression (II-11)

$$\hat{I}{}_j(\mathbf{r})\simeq \frac{2}{i\omega }{e}^{-ik{\mathring{\mathbf{\kappa }}}^i·\mathbf{r}}\sum _{\mathit{\mu }}{\hat{A}}_{\mathit{\mu }}^j(\omega ){\cos }^2\left(\frac{{\gamma }_{\mathbf{n}}^{\mp }}{2}\right){\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\mp }(\mathbf{r},\omega ),$$

(31)

where ${\gamma }_{\mathbf{n}}^{\mp }$ represents the angle between the illumination direction $-{\mathring{\mathbf{\kappa }}}^i$ and the beam direction ${\mathring{\mathbf{\kappa }}}_{\mathit{\mu }}^j$ (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\left(\mathbf{r}\right)$ and the beam data is given by (I-21)

$${\hat{A}}_{\mathit{\mu }}^j(\omega )\simeq \langle O(\mathbf{r}),{\hat{\mathsf{\Lambda }}}_{\mathit{\mu }}^j(\mathbf{r},\omega ){\rangle }_V, {\hat{\mathsf{\Lambda }}}_{\mathit{\mu }}^j(\mathbf{r},\omega )=k/2i\cos {\theta }_{\mathbf{n}}{e}^{-ik{\mathring{\mathbf{\kappa }}}^i·\mathbf{r}}{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\mp }(\mathbf{r},\omega ),$$

(32)

where the integration covers the entire scatterer domain. Thus, within the Born approximation, the data is described as projections of $O\left(\mathbf{r}\right)$ on the beam axis, using the projection kernels ${\hat{\mathsf{\Lambda }}}_{\mathit{\mu }}^j(\mathbf{r},\omega )$. As shown in I-Section 5-B in [36], this projection extracts the local stratification of O along the beam axis at an angle ${\gamma }_n^{\mp }$ defined in (31). This implies that the scattering amplitudes ${\hat{A}}_{\mathit{\mu }}^j$ are obtained from Snell’s reflections due the stratification components in $O\left(\mathbf{r}\right)$ along the $\mathit{\mu }$ beam axis, so that strong amplitudes are obtained only for those $\mathit{\mu }$ (locations and directions) that correspond to the stratification of $O\left(\mathbf{r}\right)$ along the $\mathit{\mu }$ 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 ${\hat{A}}_{\mathit{\mu }}^j$. Substituting (5) into (30) and changing the order of integrations, ${\hat{A}}_{\mathit{\mu }}^j$ can be expressed as (I-22)

$${\hat{A}}_{\mathit{\mu }}^j\simeq (\frac{k}{2\pi }{)}^2\int d^2\mathbf{\xi }\frac{ik}{2\zeta }\hat{\tilde{\phi }}{}_{\mathit{\mu }}^{\ast }(\mathbf{\xi })\overline{O}(\mathbf{K}){\mid }_{\mathbf{K}=k({\mathring{\mathbf{\kappa }}}_j-{\mathring{\mathbf{\kappa }}}^i)}.$$

(33)

The expression is the local alternative to the plane-wave-based DT identity in (8). It is recognized as ${\mathring{\mathbf{\kappa }}}_{\mathbf{n}}^{\pm }$ samples of the value of $\overline{O}\left(\mathbf{K}\right)$ over the shifted Ewald sphere defined in (9). The spectral width of these samples is that of ${\hat{\tilde{\phi }}}_{\mathit{\mu }}(\mathbf{\xi })$.

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_1$ plane, as illustrated in Figure 5a. We illustrate here only the reflected data on the $z_1$ 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_2$ is available, then the receiver array considerations are similar.

2.

The data is measured by exciting the sources one at a time by short pulse $F\left(t\right)$ that spans the desired frequency band $\Omega \in [{\omega }_{\min },{\omega }_{\max }]$ as needed to obtain the desired $K$-space coverage. The result is an $M\times M$ data matrix ${\mathrm{U}}_{pq}^s\left(t\right)$ describing the response at the p receiver due to an excitation by the q source.

3.

The data is sampled at the proper Nyquist rate and then converted to the FD via FFT, giving rise to the data matrix $\hat{\mathbf{U}}{}_{pq}^s\left(\omega \right)$. Before the calculation, the time-series are padded by zeros to avoid aliasing of the final image when it is generated by integration over all the frequencies. Note that in some applications, $\hat{\mathbf{U}}{}_{pq}^s\left(\omega \right)$ may be measure directly in the FD.

4.

The response to time-harmonic PW excitations at different angles is synthesized from $\hat{\mathbf{U}}{}_{pq}^s\left(\omega \right)$ by q-stacking the array data with proper phase terms. The result provides the PW data to the phase–space beam-based processing.

One may also calculate the time-harmonic PW spectrum of the scattered field via p-stacking with proper phase terms. The result is an $M\times M$ data matrix $\widehat{\tilde{\mathbf{U}}}{}_{p\prime q\prime }^s\left(\omega \right)$ describing the $p^{\prime }$ spectral PW due to an excitation by the $q^{\prime }$ 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.

5.

As illustrated in Figure 5a, the spectral information that can be covered by the array is determined by the size of the array and by the target range R. The size of the array should be chosen to provide sufficient spectral coverage of the target. In general, R should be within the Fresnel zone of the array, i.e., $R\ll {\left(D{\zeta }^i\right)}^2/\lambda$. Note also that due to the finite size of the array, one should avoid the array transition zone illustrated in Figure 5a.

6.

The array elements inner spacing should, in general, satisfy the Nyquist sampling rate $kd=\pi$. However, since the target range satisfies $kR\gg 1$, it is only required that the phase difference between adjacent elements will be small, yielding a sparser array with $d<\sqrt{\frac{\pi }{8k_{\max }}R}$.

Phase II—Constructing the phase space lattice

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

1.

Choosing the beam parameter b: As discussed in Section 3.2, the ID Gaussian windows in (23) are fully determined by the parameter b. The considerations of choosing b were widely studied in [21,36] for the application of local inverse scattering. This parameter balances between the beam collimation length and the beamwidth. We choose b to be on the order of the DoI domain so that the beams are collimated throughout the DoI while being small enough for transversal resolution.

2.

The phase–space lattice: The guidelines for constructing the UWB beam lattice are discussed in Section 3.2. As discussed there, we choose ${\nu }_{\max }=1/3$ which balances between stable expansion frame and moderate over-completeness (relatively small number of elements). The optimal values for $(\overline{x},\overline{\xi })$ are given by (24).

3.

The phase–space propagators: The frame elements $\left\{{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right),{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right)\right\}$ are calculated via (26). For the Gaussian windows (23), the results are ID-GB propagators (see explicit expressions in (Appendices A,C [21]).

4.

Limited physical data: We need to consider only beams whose initiation point ${\mathbf{x}}_{\mathbf{m}}$ are supported by the array size. The maximal value of ${\mathbf{\xi }}_{\mathbf{n}}$ is determined by the scan angle. If, for example, the scan angle is ${60}^{\circ }$, then $|{\mathbf{\xi }}_{\mathbf{n}}|<\frac{\sqrt{3}}{2}$.

Phase III—Calculating the beam data

1.

Calculating the expansion coefficients: The coefficients ${\hat{A}}_{\mathit{\mu }}^j$ are calculated via (30).

2.

Filtering out low amplitude data: As discussed in Section 4.2, the beam data is related to the local Snell’s reflections of the beams by the local stratification in $O\left(\mathbf{r}\right)$, which in turns are determined by the LRT of $O\left(\mathbf{r}\right)$. We therefore threshold low amplitude beams at a level of 40 dB.

Phase IV—Local reconstruction via beam backpropagation

1.

Beam backpropagation within the DoI: Next, following Section 4.1, we backpropagate the beams whose amplitudes ${\hat{A}}_{\mathit{\mu }}^j$ are larger than the threshold set above. We consider only the beams passing through the DoI (see Figure 6) or no further than 3 beam-widths away from the DoI (this distance is consistent with an effective threshold of 40 dB).

2.

The image: The imaging fields are calculated via (31), and finally the image is calculated via Equations (12)–(13).

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

The medium is plotted in Figure 7a in a $2D$ coordinate frame $\mathbf{r}=(x,z)$, with the DoI being the $20\times 20$ black rectangle. For simplicity we normalize all space-units such that the background wave-speed is $v_0=1$. Note that the contrast is relatively large with values of $O_{\max }=\pm 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 $K_z$ axis (see discussion below). Referring to the discussion in Section 2.4, it can be recovered using UWB reflection data on the $z_1$ plane. We therefore use illumination by a z-propagating time-harmonic PW with frequencies in the band $\Omega =\left[0.1,1\right]$. 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 $\left|x\right|<250$ with inter-element spacing $d=1.15$$\pi$ (larger than the Nyquist distance).

We set $b=50$, such that the beams are collimated inside the DoI, while maintaining $k_{\min }b\gg 1$ as required for collimation after (23). Using also $k_{\max }=1$ and ${\nu }_{\max }=1/3$ we obtain from (24) $(\overline{x},\overline{\xi })=(9.71,0.194)$. The resulting BF propagators $\left\{{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right),{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right)\right\}$ are calculated via (Equations (C1)–(C5) [21]).

Next we calculate the beam-domain data ${\hat{A}}_{\mathit{\mu }}^j$ via (30). The reconstructed object inside the DoI is found using the reflected data imaging field ${\hat{I}}_1\left(\omega \right)$ of (31) where we consider only backpropagated beams whose $\mathit{\mu }$ 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\pm 6$.

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 $\left|\mathbf{K}\right|\approx 0$. As discussed in Section 2.4, this missing data is due to the low frequency cutoff $k_{min}=0.1$ in the data. The main drawback of the “UWB reflection mode inversion” schemes is the missing transversal spectrum components and the $\left|\mathbf{K}\right|\to 0$ components, (which are small in this example). In general, one may try to recover this data by using transition mode data $\left(z_2\right)$ 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,43,44,50].

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., $\overline{O}\left(\right|\mathbf{K}|=0)\ne 0$). As before, we consider a $2D$ problem with $r=(x,z)$ and normalize the space-units such that the background wave-speed is $v_0=1$.

The source array is located on the $z_1=-150$ plane over $\left|x\right|<250$, with inter-element spacing $d=1.15\pi$. Using this array, we may synthesize PW illumination over a spectral range of $\pm {60}^{\circ }$ (see discussion in Section 5, Phase I(1–4)). The frequency band is $\Omega =\left[0.1,1\right]$. We consider both the reflection and transmission data over similar receiver arrays at $z_1=-150$ and $z_2=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\gg 1$ as required for collimation after (23). Using also $k_{\max }=1$ and ${\nu }_{\max }=1/3$ we obtain from (24) $(\overline{x},\overline{\xi })=(6.47,0.32)$. The resulting BF propagators $\left\{{\hat{\mathsf{\Psi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right),{\hat{\mathsf{\Phi }}}_{\mathit{\mu }}^{\pm }\left(\mathbf{r}\right)\right\}$ are calculated via (Equations (C1)–(C5) [21]).

Figure 10b depicts the reconstructed objects in the front (left) using a single PW illumination at ${\theta }_i=0$ and reflection data at the $z_1$ 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 $\left|\mathbf{K}\right|=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,43,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 ${\theta }_i=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 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 ${\theta }^i=\mp {30}^{\circ },\mp {40}^{\circ }$. The reconstruction inside the DoI is much better than in Figure 10. As before, the reconstruction outside this DoI is poor.

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:

1.

Local imaging within a given domain of interest (DoI).

2.

Reduced complexity since it accounts only for the beam basis-functions that cover the DoI.

3.

Reduced noise level since data and noise arriving from other regions are a priori filtered out.

4.

Backpropagation and imaging over a non-homogeneous background.

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.