Ultrawideband Microwave Sensing and Imaging Using Time-Reversal Techniques: A Review

Abstract

This paper provides an overview of some time-reversal (TR) techniques for remote sensing and imaging using ultrawideband (UWB) electromagnetic signals in the microwave and millimeter wave range. The TR techniques explore the TR invariance of the wave equation in lossless and stationary media. They provide superresolution and statistical stability, and are therefore quite useful for a number of remote sensing applications. We first discuss the TR concept through a prototypal TR experiment with a discrete scatterer embedded in continuous random media. We then discuss a series of TR-based imaging algorithms employing UWB signals: DORT, space-frequency (SF) imaging and TR-MUSIC. Finally, we consider a dispersion/loss compensation approach for TR applications in dispersive/lossy media, where TR invariance is broken.

1. Introduction

Remote sensing systems in the microwave and millimeter wave (mm-wave) frequency ranges provide unique capabilities for detection and imaging of obscured targets. For example, microwaves can penetrate through foliage to detect obscured vehicles and personnel in forest environments. Mm-waves can penetrate through non-conductive walls and many packaging. They can also penetrate through dust, fog and smoke (as opposed to optical- and infrared-based systems). Even for scenarios where there is no distinctive advantage among various sensor modalities, microwave and mm-wave systems can play an important role because of their ease of integration into multimodal sensing apparatus [1].

Challenges for the full exploitation of microwave and mm-wave frequencies include vulnerability to certain atmospheric and meteorological phenomena and attenuation inside lossy materials. In addition, intervening media are often disordered with complex constitutive properties and/or many secondary scatterers about which precise information is not available. As a result, signals from the primary target(s) are often weak and/or masked by clutter and multipath (multiple scattering) contributions, which confounds detection, causes erratic tracking and makes it difficult to extract relevant information using conventional imaging and classification algorithms.

However, recent results have shown that multipath can actually be exploited to improve detection and imaging capabilities in sensing scenarios (or, somewhat equivalently, channel capacity in wireless communication scenarios) [2,3]. One such technique is time-reversal (TR) [4,5,6,7], which exploits multipath components in the intervening media to achieve superresolution, i.e., resolution that beats the classical diffraction limit. TR originated in acoustics and relies on the TR invariance of the wave equation in stationary and lossless media. TR involves physical or synthetic retransmission of signals acquired by a set of transceivers in a time-reversed fashion, i.e., last-in first-out. The retransmitted signals propagate “backwards”, naturally reversing the path that they underwent during forward propagation, which results in (automatic) energy focusing around initial source location(s). The “source” in this case can be either active (transmit mode) or passive (scatterer or “echo” mode).

The TR techniques that rely on ultrawideband (UWB) operation are further attractive because they can exploit advantages of simultaneous operation at low (e.g., more penetration into lossy materials) and high frequencies (e.g., better resolution), and because they enable imaging techniques in random media that depend only on the statistical properties (instead of a particular realization) of the random media, i.e., they are statistically stable [8,9,10].

The first successful TR experiments were carried by Fink and his collaborators using acoustic waves [11]. The underlying properties of TR operation were later investigated in more details both theoretically [10,12,13,14,15] and numerically [8,9,16,17,18,19].

Since then, TR techniques have experienced explosive growth in applications. Among the TR applications in medicine, one can cite lithotripsy [20,21], ultrasonic focusing through the skull [22], hyperthermia [23,24], remote inspection of internal masses [25] and health monitoring [26,27]. Applications of TR in geophysics and geoscience have included probing earthquake epicenters [28] and detecting objects buried underground [29]. Underwater applications have included sonar and acoustic communications [30,31], intruder detection, and echo-to-reverberation enhancement [32,33]. A commercial phased array sound system utilizing time-reversal concept is also developed in [34].

Nondestructive testing applications have utilized it for detecting defects in materials and structures [35,36,37,38,39].

In addition, imaging using synthetic TR techniques [40,41], TR interferometry [42,43,44] and imaging for sensor networks [45] have been addressed for both acoustic and electromagnetic waves. For most TR applications involving imaging, analysis of the so-called TR operator (TRO) is fundamental. In particular, the eigenspace analysis of the TRO provides important information about the scattering scenario under study. In the case of multiple targets for example, selective focusing using successive TR waves becomes possible. Perhaps the most basic approach that exploits the eigenvalue/vector structure of the TRO is the “DORT” (a French language acronym standing for “decomposition of the time-reversal operator") method [29,46,47,48,49,50,51]. In a similar vein, the “TR-MUSIC” method (derived from the MUltiple SIgnal Classification algorithm widely used for angle-of-arrival estimation in radar applications) has also been developed [40,52,53,54] and considered for a variety of applications [53,54,55,56,57].

In the frequency-domain, TR is equivalent to a phase-conjugation which has found applications in optics [58,59,60] as well as in electromagnetics where it is used to cancel distortions in the medium and provide beam steering action by using self-phasing or retrodirective array antennas [61,62,63,64,65]. In [66], it was experimentally shown that by using a single frequency phase-conjugation approach, energy focusing around two targets is also possible.

Among the most notable TR experiments using RF and microwaves, one can cite [67,68] where it was experimentally shown that a 1 MHz wideband pulse can be focused in a cavity environment.

Various TR-based target/change detection and interference cancelation algorithms have been discussed and investigated experimentally in [69,70,71,72,73].

Breast cancer detection using microwaves has been a topic of much interest [74]. TR is quite suited for this problem because of its natural ability to focus energy on malignant tissues (stronger scatterer). TR-based approaches for detection of breast cancer using microwaves are discussed in, e.g., [75,76,77,78,79,80,81].

More generally, TR-based electromagnetic imaging and target detection in cluttered environments with discrete scatterers have been experimentally investigated in [82,83,84].

The DORT method has been discussed in connection with remote sensing of buried objects in [29,49,85]. Additionally, the TR operator for electromagnetic waves in homogeneous medium was recently analyzed in [86,87].

TR-enabled focusing in a forest environment has been investigated numerically in [88,89]. Moreover, a TR synthetic aperture radar (TR-SAR) imaging method was developed and applied to an environment with strong multipath in [90,91]. Recently filed patent applications consider TR-based techniques for SAR images with improved target focusing and ghost image removal [92] as well as beamforming [93]. Additionally, TR technique has been applied for the through wall imaging in [94].

Finally, TR techniques have also found a number of applications in wireless communications such as the development of TR-based spatio-temporal matched filters to reduce channel dispersion and inter-symbol interference, thereby increasing the channel capacity [70,95,96,97,98,99,100,101,102]. Finally, a practical implementation of time-reversal of broadband microwave signals using optics has been demonstrated in [103] as well as in [104].

The objective of this review paper is to provide a summary of some TR-based techniques for detection and imaging applications using UWB microwaves. We first discuss TR operation in continuous random media which acts as a scenario where superresolution effects are manifested [105]. Then, we examine two examples from the “signal space” family of TR methods, viz. the time-domain (TD) DORT method [106], which allows selective focusing of multiple scatterers, and the space-frequency UWB TR imaging method [107], which utilizes spatial and frequency information of the received signals in tandem. Next, we discuss the “null space” TR-MUSIC method [108], which outperforms signal space methods under weak clutter conditions. Finally, we illustrate a dispersion/loss compensation approach aimed at extending TR operation to dispersive and lossy media, where TR invariance is broken [109].

2. Time-Reversal of Electromagnetic Waves and Superresolution

The TR concept exploits the invariance of the wave equation under TR, i.e., if $\overline{E}(\overline{r},t)$ is a solution of the wave equation in lossless media,

$${\nabla }^2\overline{E}(\overline{r},t)-\mu \left(\overline{r}\right)ϵ\left(\overline{r}\right)\frac{{\partial }^2}{\partial t^2}\overline{E}(\overline{r},t)=0$$

(1)

then $\overline{E}(\overline{r},-t)$ is also a solution. Here, $\overline{E}(\overline{r},t)$ represents the electric field, $\overline{r}=\widehat{x}x+\widehat{y}y+\widehat{z}z$ is the position vector, and μ and ϵ are the permeability and permittivity of the medium, respectively. This property of the wave equation guarantees that for every wave propagating away from the source, there exists a reversed wave that would precisely retrace the path of the original wave back to the source. This occurs regardless of the presence of scattering objects and of inhomogeneities in μ and ϵ. The reversed waves would converge and focus coherently at the original source location as if time were flowing backwards, $t\to -t$.

This concept has been exploited by recording the wave using a set of transceivers and sending it back to its source as if time had been reversed, thereby yielding a variety of applications as mentioned in the Introduction. In practice, it is often not possible to completely surround the source so that TR is more generally performed with a set of receive/transmit antennas, the so-called TR array (TRA), that comprises only a limited angle. This causes some “information loss”; in general, the focal spot size gets larger as the array aperture size becomes smaller. Note that even though the term “array” is predominantly used in the TR literature, the antenna elements in the TRA ideally should work as independent transceivers, more akin to a “MIMO (multiple-input multiple-output) radar”. Recovering some of the “lost information” is possible through the exploitation of multipath components, as illustrated in Figure 1.

Figure 1. Time reversal experiment using a limited aspect array. (left) Forward Propagation of the short input pulse, (right) Backward Propagation of the time-reversed signals.

Figure 1. Time reversal experiment using a limited aspect array. (left) Forward Propagation of the short input pulse, (right) Backward Propagation of the time-reversed signals.

Here, a point source located at ${\overline{r}}_0$ transmits a short UWB pulse $s\left(t\right)$. The transmitted signal propagates through the medium and is received by an antenna array. The signal received at the $i^{th}$ antenna is

$$f_i\left(t\right)=s\left(t\right){*}_t{h}_{{\overline{r}}_0{\overline{r}}_i}\left(t\right)$$

(2)

where ${*}_t$ denotes time convolution and $h_{{\overline{r}}_0{\overline{r}}_i}\left(t\right)$ is the impulse response between the antennas at ${\overline{r}}_0$ and ${\overline{r}}_i$. Reciprocity allows us to write $h_{{\overline{r}}_0{\overline{r}}_i}\left(t\right)=h_{{\overline{r}}_i{\overline{r}}_0}\left(t\right)$. The signals received at each array element are recorded, reversed in time, and transmitted back to the same medium (Figure 1). The time-reversed signal at the original source point ${\overline{r}}_0$ due to $i^{th}$ antenna is then given by

$$p_i\left(t\right)=\underset{f_i(-t)}{\underbrace{s(-t){*}_t{h}_{{\overline{r}}_0{\overline{r}}_i}(-t)}}{*}_t{h}_{{\overline{r}}_i{\overline{r}}_0}\left(t\right)$$

(3)

where the last two terms ($h_{{\overline{r}}_0{\overline{r}}_i}(-t){*}_t{h}_{{\overline{r}}_i{\overline{r}}_0}\left(t\right)$) represent a correlation filter (time-correlator). This correlation function has a maximum at $t=0$ which corresponds to the energy of $h_{{\overline{r}}_0{\overline{r}}_i}\left(t\right)$, i.e., $\int |h_{{\overline{r}}_0{\overline{r}}_i}{\left(t\right)|}^{2}dt$. With multiple antennas, TR system performance improves since each antenna will have a maximum at the original source location and they will constructively interfere to improve the TR peak signal. This coherent interference does not occur arbitrarily, but always at the original source location. For a TRA with N elements, the received signal at the original source location becomes:

$$p({\overline{r}}_0,t)=\sum _{i=1}^{N}s(-t){*}_t{h}_{{\overline{r}}_0{\overline{r}}_i}(-t){*}_t{h}_{{\overline{r}}_i{\overline{r}}_0}\left(t\right)$$

(4)

In addition to being a time-correlator, TR also acts as a space-correlator. In the above analysis, the TR waveform is exactly matched to the original source point ${\overline{r}}_0$. However, at any other point $\overline{r}$ in the domain, this signal can be written as

$$p(\overline{r},t)=\sum _{i=1}^{N}s(-t){*}_t{h}_{{\overline{r}}_0{\overline{r}}_i}(-t){*}_t{h}_{{\overline{r}}_i\overline{r}}\left(t\right)$$

(5)

As the probe antenna location $\overline{r}$ gets further away from the original source location ${\overline{r}}_0$, then, similar to time-correlation analysis, uncorrelated terms tend to cancel each other. For media with rich multipath components, correlation peak gets sharper and a better (sharper) focusing both in time and space can be achieved.

It is well known that the focusing spot size is dictated by the classical diffraction limit [110] which states that in an homogeneous media, cross-range ($d_c$) resolution is given by

$$d_c=\lambda \frac{L}{a}$$

(6)

where λ is the wavelength, a is the antenna aperture length and L is the distance between the array and the source (Figure 1). In scenarios with multipath, a should be replaced by the effective aperture length $a_e$ as the inhomogeneities in the intervening medium affect the receiving pattern of the TRA. As long as some of the diverging wave components are redirected toward the TRA, the effective aperture length increases ($a_e>a$) resulting in a better focusing resolution than the homogeneous medium case (superresolution) [4]. This is as illustrated in Figure 2. We will next focus on superresolution effects created by the volumetric scattering in a continuous random media. The following subsection introduces the random medium model used in this work.

Figure 2. Effective aperture increase in media with multipaths. (a) Homogeneous medium with no multipaths ($a_e=a$), (b) multipaths created by the waveguide-like structure composed of two perfect electric conductor (PEC) walls ($a_e>a$), (c) multipaths created by discrete scatterers ($a_e>a$).

Figure 2. Effective aperture increase in media with multipaths. (a) Homogeneous medium with no multipaths ($a_e=a$), (b) multipaths created by the waveguide-like structure composed of two perfect electric conductor (PEC) walls ($a_e>a$), (c) multipaths created by discrete scatterers ($a_e>a$).

2.1. Physics-Based Ultrawideband Clutter Models

Complex natural media such as snow, vegetation, rocks, soils, and some biological tissues often cannot be described in a deterministic manner. Therefore, statistical models (random medium models) should be employed instead [111]. A random or disordered media can be classified either as (i) discrete random media, characterized as a discrete set of scatterers (e.g., trees, obstacles, buildings) at random locations, or as (ii) continuous random media, characterized by pointwise fluctuations on its properties (e.g., some biological media, soils, smoke) described in terms of a stochastic process with given first- and second-order statistics (spatial correlation functions).

With subsurface sensing scenarios in mind, we focus our attention on the characterization of continuous random media using constitutive parameters from particular soil models. Soil presents natural variability in density, composition and moisture that affect its permittivity and conductivity. Additionally, since in a host of applications the aim is to detect man-made buried objects, the soil has usually been excavated and therefore it is not expected to have a homogeneous or even layered distribution. In the absence of experimental data to support a specific choice of random medium model, continuous random medium models with Gaussian distributions are preferred for their generality and mathematical properties requiring few statistical parameters. Such random media are also characterized by a spatial correlation function. The relative permittivity of the medium is described as $ϵ\left(\overline{r}\right)={ϵ}_m+{ϵ}_f\left(\overline{r}\right)$ where ${ϵ}_m$ is the mean value of the relative permittivity and ${ϵ}_f\left(\overline{r}\right)$ is a function of position characterizing the random fluctuation on $ϵ\left(\overline{r}\right)$ with $〈{ϵ}_f\left(\overline{r}\right)〉=0$. At every point in space, the fluctuation term is a Gaussian random variable with zero mean and probability density function (pdf) given by:

$$P_{{ϵ}_f}\left(\zeta \right)=\frac{1}{\sqrt{2\pi \delta }}exp\left(\frac{-{\zeta }^2}{2\delta }\right)$$

(7)

The underlying assumption with this pdf is that the fluctuating dielectric permittivity is real. However, for more general cases where the dielectric permittivity is complex and random, a similar pdf can independently be used for both the real and imaginary parts of the dielectric permittivity. But, in this paper, we only consider lossless random medium and we refer the interested readers to another paper where we have considered both random and lossy (dispersive) media [106]. The random medium is characterized by transverse and vertical correlation lengths ($l_s$ and $l_z$) and variance (δ), with the correlation function between the permittivities at two points also described by a Gaussian function

$$C({\overline{r}}_1-{\overline{r}}_2)=\langle {ϵ}_f\left({\overline{r}}_1\right){ϵ}_f^{*}\left({\overline{r}}_2\right)\rangle =\delta exp\left(-\frac{|{\overline{x}}_1-{\overline{x}}_2{|}^2+{|{\overline{y}}_1-{\overline{y}}_2|}^2}{l_s^2}-\frac{|z_1-z_2{|}^2}{l_z^2}\right)$$

(8)

Previously, similar correlation functions have been used in [111,112,113]. Throughout this paper, we assume $l_s=l_z$ so that we can isolate the effect of single correlation length on the TR focusing. Following these definitions, the procedure discussed in [114] can be used to generate the random medium realizations for numerical simulations. For modeling the scattering of embedded discrete target(s) in such media, we employ the finite-difference time-domain (FDTD) algorithm to solve Maxwell equations. The FDTD computational grid is truncated by perfectly matched layers (PML) [115] via stretched coordinates [116] to provide the necessary absorption of outgoing waves reaching the computational boundaries. A FDTD computational domain with $N_x\times N_y=200\times 240$ grid points and uniform space discretization size ${\Delta }_x={\Delta }_y={\Delta }_s=0.0137$ m is utilized in the examples considered below. A linear TRA of $N=7$ uniformly spaced dipole antennas are located just above a lossless random medium with spatially fluctuating permittivity. The TRA lies parallel to the x-direction and the dipoles are separated by ${\lambda }_c/2$ where ${\lambda }_c$ corresponds to the wavelength at the central frequency for the mean ground permittivity. The location of the central antenna is assumed as origin, i.e., ${\overline{R}}_4=(0,0){\Delta }_s$ where ${\overline{R}}_i$ is the location of the $i^{th}$ TRA antenna. An electric dipole located at ${\overline{r}}_s=(x_s,y_s)=(0,155){\Delta }_s=(0,2.12)m$ is initially fed by the current source $\overline{J}({\overline{r}}_s,t)=\widehat{e}{s}_{}\left(t\right){\delta }_d\left({\overline{r}}_s\right)$ where $\widehat{e}$ ($\widehat{x}$, $\widehat{y}$ or $\widehat{z}$) is the unit vector representing the dipole polarization, ${\delta }_d\left(\overline{r}\right)$ is the Dirac delta function and $s_{}\left(t\right)$ is the UWB time-domain excitation taken as the first derivative of the Blackmann–Harris (BH) pulse [117] that vanishes after a time period of $T=1.55/f_c$, where $f_c=400$ MHz is the central frequency.

2.2. Numerical Results

In this section, we illustrate the superresolution effects versus random medium statistical parameters, viz. variance (δ) and correlation length ($l_s$).

We first investigate the effect of the variance of the random medium on the refocusing properties of the time reversed signals for the $T{M}_z$ polarization case. The variance (δ) is varied from $0.025{ϵ}_m$ to $0.125{ϵ}_m$, while the correlation length is kept fixed at $l_s=8{\Delta }_s$. The snapshots of the z-component of the electric field ($E_z$) distribution at the time of refocusing are shown in Figure 3.

It is observed that as the variance increases, the amplitude of the focused field increases and the spot size of the focused field is reduced (somewhat counter-intuitively), characterizing superresolution. This can be explained from the fact that as the variance increases more multipath is produced. For the time reversed signals, the multipath tends to interfere constructively only at the focusing point where coherent perfect phase conjugation exists. This can also be understood as an increase on the effective aperture of the TRA due to multipathing [4,105]. Figure 5(a) shows the spatial distribution of the field components of Figure 3 at the time of refocusing at the source plane ($y=y_s$) with respect to the (transverse) x-coordinate (cross-range).

Figure 3. Spatial distribution of the time-reversed $E_z$ field component of the electric field at the time of refocusing for increasing δ and fixed $l_s$ of $l_s=8{\Delta }_s$. Plots are given in linear scale.

Figure 3. Spatial distribution of the time-reversed $E_z$ field component of the electric field at the time of refocusing for increasing δ and fixed $l_s$ of $l_s=8{\Delta }_s$. Plots are given in linear scale.

Figure 4. Cross-range resolution for varying (a) 1st and (b) 2nd order medium statistics.

Figure 4. Cross-range resolution for varying (a) 1st and (b) 2nd order medium statistics.

Figure 5. (a) Comparison of incoherent and coherent time-domain signals. Note that although (b) the magnitudes of the frequency domain representation of both signals are same, (c) the phases are not. Incoherency comes from the oscillations in the phase arising from the SVD or EVD algorithm.

Figure 5. (a) Comparison of incoherent and coherent time-domain signals. Note that although (b) the magnitudes of the frequency domain representation of both signals are same, (c) the phases are not. Incoherency comes from the oscillations in the phase arising from the SVD or EVD algorithm.

Similar simulations are performed to assess the effect of varying the correlation length on the resolution for the $T{M}_z$ polarization. The correlation length is varied from $l_s=8{\Delta }_s$ to $l_s=20{\Delta }_s$, for fixed variance of $\delta =0.04{ϵ}_m$. Cross-range resolutions are shown in Figure 5(b), where it is seen that increasing the correlation length degrades the focusing properties. This can be explained from the fact that, for larger correlation lengths, the random medium starts to behave (locally) closer to an homogeneous medium and volumetric scattering effects are reduced. Note that this observation applies only for the range of correlation values and problem size considered here. If the correlation length becomes very small versus the wavelength, then a homogenization would be applicable. In general, the effect of changes on the correlation length on the focusing resolution is more pronounced when the correlation length is on the order of the wavelength. Another point to consider is the effect of polarization. As shown in [105], increased multiple scattering tends to increase depolarization, hence fully polarimetric data should be utilized for improved refocusing resolution. We also note that, although only the results of single realizations are shown here, another interesting phenomenon that occurs in connection with time-reversed back-propagation of wavefields in a random medium is the occurrence of the so-called time-domain statistical stability [8,9]. This means that the retrofocused field is self-averaging in the time-domain for a single realization and does not depend on the particular realization of the random medium. Self-averaging occurs because the time reversed field is equivalent to a phase conjugated field. After back-propagation, the conjugate phase exactly cancels out the random phase of the initial field and different frequency components add fully coherently only nearby the original source location, hence the focal spot is not affected by the particular realization.

3. Time-Reversal based Detection and Imaging Methods

In this section, we review some TR-based signal processing techniques to achieve UWB EM detection and imaging for multiple objects (targets) embedded in inhomogeneous media. This requires the analysis of the eigenspace of the TR operator (TRO) and utilization of its signal or null spaces. The TRO is obtained from the multistatic data matrix (MDM). The MDM element ${\left[\mathit{K}\left(\omega \right)\right]}_{ij}$ corresponds to the signal received at the ith antenna when the jth antenna is transmitting.

3.1. Signal Space Methods

3.1.1. Time-Domain DORT

The DORT method considers a $N\times N$ MDM $\mathit{K}\left(\omega \right)$ (where ω is the frequency) obtained by probing the medium by a TRA of N transceivers. Since in this case each MDM element is associated with a different spatial location, this type of MDM is labeled as “space-space MDM”. In the frequency domain, TR of $\mathit{K}\left(\omega \right)$ is represented by its Hermitian conjugate ${\mathit{K}}^{†}\left(\omega \right)$ and the TRO is defined by $\mathit{T}\left(\omega \right)={\mathit{K}}^{†}\left(\omega \right)\mathit{K}\left(\omega \right)$. The SVD of the MDM is given by

$$\mathit{K}\left(\omega \right)=\mathit{U}\left(\omega \right)\Lambda \left(\omega \right){\mathit{V}}^{†}\left(\omega \right)$$

(9)

where $\mathit{U}\left(\omega \right)$ and $\mathit{V}\left(\omega \right)$ are unitary matrices and $\Lambda \left(\omega \right)$ is real diagonal matrix of singular values. The EVD of the TRO can be written as $\mathit{T}\left(\omega \right)=\mathit{V}\left(\omega \right)\mathit{S}\left(\omega \right){\mathit{V}}^{†}\left(\omega \right)$, where $\mathit{S}\left(\omega \right)={\Lambda }^{†}\left(\omega \right)\Lambda \left(\omega \right)$ is the diagonal matrix of eigenvalues. The columns of the unitary matrix $\mathit{V}\left(\omega \right)$ are normalized eigenvectors of the TRO (${\mathbf{v}}_i\left(\omega \right),i=1,..,N$). The DORT method utilizes the eigenvectors of the TRO signal space (SS) which is formed by the eigenvectors with non-zero (significant) eigenvalues, i.e.,

$$SS\left(\omega \right)=\{{\mathbf{v}}_1\left(\omega \right),...,{\mathbf{v}}_{M_s}\left(\omega \right)\}\mathrm{with}{\lambda }_1^2\left(\omega \right)>..>{\lambda }_{M_s}^2\left(\omega \right)>0$$

(10)

where $M_s$ is the number of significant eigenvalues. For isotropic scattering from well-resolved point-like scatterers, each significant eigenvalue of the TRO is associated with a particular scatterer. Subsequent back-propagation of the corresponding eigenvector yields a wavefront focusing on that particular scatterer only [118]. Therefore, selective focusing on the mth scatterer is achieved by exciting the TRA with $N\times 1$ column vector ${\mathbf{r}}_m\left(\omega \right)$ generated from ${\lambda }_m\left(\omega \right)$ and eigenvector ${\mathbf{v}}_m\left(\omega \right)$ via

$${\mathbf{r}}_m\left(\omega \right)={\mathit{K}}^{†}\left(\omega \right){\mathbf{u}}_m\left(\omega \right)={\lambda }_m\left(\omega \right){\mathbf{v}}_m\left(\omega \right)$$

(11)

This represents the single-frequency DORT method. For UWB signals, EVD can be applied at all the available frequencies and a time-domain signal can be generated by

$${\mathbf{r}}_p\left(t\right)={\mathcal{F}}^{-1}\left({\mathit{K}}^{†}\left(\omega \right){\mathbf{u}}_p\left(\omega \right)\right)={\mathcal{F}}^{-1}\left({\lambda }_p\left(\omega \right){\mathbf{v}}_p\left(\omega \right)\right)$$

(12)

where ${\mathcal{F}}^{-1}$ denotes the inverse Fourier transformation. Backpropagation of these time-domain signals into the probed medium characterizes the TD-DORT method [106]. In cases where the background Green’s function is not known (which is the case for most subsurface sensing scenarios), an approximate Green’s function ($G_0$) can be used to obtain synthetic images of the probed medium by using the following TD-DORT imaging functional:

$$\begin{array}{ccc}\hfill {\mathbf{D}}_p^{\Omega }\left({\overline{X}}_s\right)& =& 〈{\mathbf{g}}_0({\overline{X}}_s,t),{\mathbf{r}}_p\left(t\right)〉{\mid }_{t=0}\hfill \\ & =& \sum _{i=1}^N{r}_p^i(-t){*}_t{G}_0({\overline{X}}_s,{\overline{R}}_i,t){\mid }_{t=0}\hfill \\ & =& {\int }_{\Omega }{\lambda }_p\left(\omega \right){\mathbf{v}}_p^T\left(\omega \right){\mathbf{g}}_0({\overline{X}}_s,\omega )d\omega \hfill \end{array}$$

(13)

where

$${\mathbf{g}}_0({\overline{X}}_s,t)={\left[G_0({\overline{X}}_s,{\overline{R}}_1,t),\dots ,G_0({\overline{X}}_s,{\overline{R}}_N,t)\right]}^T$$

(14)

is the approximate time-domain steering vector connecting any search point ${\overline{X}}_s$ in the probed domain to the antenna locations ${\overline{R}}_i$ for $i=1,..,N$, $r_p^i$ and $v_p^i$ are the ith element of ${\mathbf{r}}_p$ and ${\mathbf{v}}_p$, respectively and Ω is the bandwidth of operation. If one assumes, for example, that at least the average constitutive parameters of the random medium are known with a good degree of confidence, then $G_0$ is chosen the Green’s function of an homogeneous medium with the mean constitutive parameters of the underlying random (inhomogeneous) medium that produces $\mathit{K}\left(\omega \right)$. Note also that when the background medium is not known, energy detector which is a non-coherent imaging method can also be utilized [72].

For well-resolved scatterers, the cross-range resolution of the DORT method is related to the classical diffraction limit (Equation 6) which is directly proportional to the wavelength and propagation distance, and inversely proportional to the effective aperture length of the TR array. When the well-resolvedness criterion is broken, eigenvectors of the signal space of space-space MDMs might become linear combinations of the Green’s function vectors connecting the targets to the TR array, which results in a degradation (blurring) of the image. Another problem with the time-domain DORT stems from the eigenvalue decomposition step which creates eigenvectors with arbitrary frequency dependent phase ${\varphi }_{svd}\left(\omega \right)$. Direct combination of these eigenvectors creates incoherent time-domain signals which in turn severely affect the TD-DORT excitation signals. This is particularly important in inhomogeneous background media with strong multiple scattering where multipath components need to be coherently combined over the entire bandwidth. Therefore, a pre-processing step should be applied to the space-space eigenvectors to obtain coherent time domain signals by canceling the arbitrary phase term ${\varphi }_{svd}\left(\omega \right)$. Such pre-processing steps can be done by projecting the incoherent eigenvectors onto the columns of the MDM [8] or by a phase smoothing algorithm which tracks the phase difference between the adjacent frequency eigenvectors [106,119,120]. The time- and frequency-domain signals obtained for a random medium case with and without phase smoothing algorithm are shown in Figure 6(a). The effectiveness of the phase smoothing algorithm is evident from these plots.

Figure 6. The first eigenvalue distribution (most and only significant one in this case) of the space-space MDM with respect to the frequency (left) and singular values of the individual (center) and full (right) space-frequency MDMs.

Figure 6. The first eigenvalue distribution (most and only significant one in this case) of the space-space MDM with respect to the frequency (left) and singular values of the individual (center) and full (right) space-frequency MDMs.

3.1.2. Space-Frequency TR Imaging

An alternative approach to pre-processing is based on the application of the SVD to a different type of MDM which incorporates both sensor location and UWB frequency data simultaneously. Such MDM is denoted as “space-frequency (SF)” MDM. This approach was first introduced in connection with MUSIC-type algorithms [121]. Here, we illustrate them in the context of TR-based UWB imaging algorithms.

SF-MDMs can be obtained by transmitting an UWB pulse from the $n^{th}$ antenna and recording the received signals from all the TRA receiver antennas to yield an $N\times M_f$ matrix given by

$${\mathit{K}}_{SF}^n=\left(\begin{array}{ccc}{k}_{1n}\left({\omega }_1\right)& \cdots & k_{1n}\left({\omega }_{M_f}\right)\\ ⋮& \cdots & ⋮\\ k_{Nn}\left({\omega }_1\right)& \cdots & k_{Nn}\left({\omega }_{M_f}\right)\end{array}\right)$$

(15)

where each row consists of the uniform samples of the Fourier transform of the time-domain signal corresponding to the respective MDM element, and $M_f$ is the number of frequency samples. Once N of these individual SF-MDMs are obtained, SVD is applied to each of them to yield ${\mathit{K}}_{SF}^n={\mathit{U}}_{SF}^n{\Lambda }_{SF}^n{\left({\mathit{V}}_{SF}^n\right)}^{†}$, where ${\mathit{U}}_{SF}^n$ is the $N\times N$ matrix of left singular vectors, ${\mathit{V}}_{SF}^n$ is the $M_f\times M_f$ matrix of right singular vectors, and ${\Lambda }_{SF}^n$ is the $N\times M_f$ matrix of singular values. Note that ${\mathit{K}}_{SF}^n$ maps the frequency space onto the receiver space ${\mathit{K}}_{SF}^n{\mathbf{v}}_{S{F}_i}^n={\lambda }_{S{F}_i}^n{\mathbf{u}}_{S{F}_i}^n$, where ${\lambda }_{S{F}_i}^n$ is the $i^{th}$ singular value, ${\mathbf{v}}_{S{F}_i}^n$ is the $i^{th}$$M_f\times 1$ right singular vector that represents the frequency content of the received signals, and ${\mathbf{u}}_{S{F}_i}^n$ is the $i^{th}$$N\times 1$ left singular vector containing spatial (sensor location) information. The left singular vectors ${\mathbf{u}}_{S{F}_i}^n$ for $i=1,\dots ,N$ form an orthonormal set spanning the sensor location space; similarly, right singular vectors ${\mathbf{v}}_{S{F}_i}^n$ for $i=1,\dots ,M_f$ form an orthonormal set spanning the frequency space. Inverse Fourier transformation can be applied to the right singular vectors to obtain time-domain signals which are coherent. In other words, SVD applied to the SF-MDM does not create arbitrary phase dependent term as seen in the TD-DORT implementation. Therefore, the time-domain excitation signals to be backpropagated can be approximated as

$$s^n\left(t\right)=\sum _{k=1}^P{\lambda }_{S{F}_k}^n{v}_{S{F}_k}^n\left(t\right)$$

(16)

where $v_{S{F}_k}^n\left(t\right)={\mathcal{F}}^{-1}\left(v_{S{F}_k}^n\left(\omega \right)\right)$ and P is the total number of time-domain signals being included in the approximation. This is determined by examining the singular values and the associated singular vectors. Time-domain signals corresponding to relatively small singular values and those with erratic behavior are not included. The resulting time-domain signals provide the UWB frequency data but they do not possess any sensor location data. Therefore, we use the left singular vectors to provide the necessary amplitude and phase shifts to be applied to each TR antenna during backpropagation. To this end, we define a vector functional

$$\mathbf{f}(\mathbf{a},z\left(t\right))={\mathcal{F}}^{-1}\left\{{[A_0{e}^{j{\varphi }_0}z\left(\omega \right),\dots ,A_N{e}^{j{\varphi }_N}z\left(\omega \right)]}^T\right\}$$

(17)

where $z\left(t\right)$ is the time-domain signal to be used and $\mathbf{a}={[A_0{e}^{j{\varphi }_0},\dots ,A_N{e}^{j{\varphi }_N}]}^T$ is the $N\times 1$ vector to determine relative time-delays and amplitudes. This functional yields a time-domain vector ${\mathbf{r}}_{S{F}_i}^n\left(t\right)=\mathbf{f}\left({\mathbf{u}}_{S{F}_i}^n,s^n\left(t\right)\right)$ which is to be backpropagated from the receivers. If the background medium is not known precisely, an approximate Green’s functions can be utilized, which results in the following imaging functional:

$$\begin{array}{ccc}\hfill {\mathbf{I}}_{{SF}_i}^n\left(\overline{X_s}\right)& =& \langle {\mathbf{g}}_0({\overline{X}}_s,t),{\mathbf{r}}_{S{F}_i}^n\left(t\right)\rangle {\mid }_{t=0}\hfill \\ & =& \sum _{k=1}^N{r}_{S{F}_{ik}}^n(-t){*}_t{G}_0({\overline{X}}_s,{\overline{R}}_k,t){\mid }_{t=0}\hfill \\ & =& {\int }_{\Omega }{\left({\mathbf{r}}_{S{F}_i}^n\left(\omega \right)\right)}^{†}{\mathbf{g}}_0({\overline{X}}_s,\omega )d\omega \hfill \\ & =& {\int }_{\Omega }{\left(s^n\left(\omega \right)\right)}^{*}{\left({\mathbf{u}}_{S{F}_i}^n\right)}^{†}{\mathbf{g}}_0({\overline{X}}_s,\omega )d\omega \hfill \end{array}$$

(18)

where * represents complex conjugation. This procedure is repeated for all $n=1,\dots ,N$. The resulting image is averaged via

$${\mathbf{I}}_{S{F}_i}\left(\overline{X_s}\right)=\frac{1}{N}\sum _{n=1}^N{\mathbf{I}}_{{SF}_i}^n\left(\overline{X_s}\right)$$

(19)

and denoted as “SF-image”. Note that instead of the left singular vectors, one can also utilize the eigenvectors obtained from the DORT method at the central frequency for the spatial information [107]. A second kind of SF-MDM can be obtained by combining all the ${\mathit{K}}_{SF}^n$ for $n=1,..,N$ into a single $N^2\times M_f$ matrix given as

$${\mathit{K}}_{SF}^{full}={[{\mathit{K}}_{SF}^1;{\mathit{K}}_{SF}^2;\dots ;{\mathit{K}}_{SF}^N]}^T$$

(20)

This is denoted as the “full SF-MDM”. Similar to before, SVD applied to this matrix yields $N^2\times N^2$ left and $M_f\times M_f$ right singular vectors. The right singular vectors of the full SF-MDM are very similar to the individual SF-MDM and can be utilized to obtain time-domain excitation signals for backpropagation. The sub-vectors of the left singular vectors can provide the sensor location data for the SF imaging. For further details, readers are referred to [107]. Finally, we note that the phase and magnitude of the significant left singular vectors highly depend on the spatial distribution of the TR array antennas (sensor). As a result, inter-distances between sensors affect the phase and magnitude distributions of the left singular vectors. Similarly, for the DORT method, it will affect the phase and magnitudes of the eigenvectors. Depending on the scenario, this change can have positive or negative effects on the final imaging performance and has to be studied in more details which is reserved as a future study.

Next, we apply both TD-DORT and SF-imaging to a subsurface detection scenario to illustrate the similarities and differences between the two methods.

3.1.3. Results

To demonstrate the performance of DORT and SF-imaging methods, we employ the same simulation scenario as considered before except for the following differences: (a) we assume the central antenna located at the origin (i.e. $R_4=(0,0){\Delta }_s$) and (b) a single scatterer at $(30,80){\Delta }_s$ is considered in both homogeneous and random media, the latter characterized by variance $\delta =0.03$ and correlation length $ls=6{\Delta }_s$. Figure 6 shows the corresponding eigenvalues for the TD-DORT method and singular values for the SF-imaging method.

It is observed that with increasing multiple scattering (variance), the magnitude of the eigenvalues and singular values increases. Additionally, the first dominant singular values obtained in homogeneous (HM) and random medium (RM) are close to each other and exhibit similar time-domain signatures (except for fluctuations due to clutter in the random medium case) as seen in Figure 7.

Figure 7. The first two significant time-domain right singular vectors obtained in homogeneous and random media using both the individual and full SF-MDM.

Figure 7. The first two significant time-domain right singular vectors obtained in homogeneous and random media using both the individual and full SF-MDM.

Among the remaining non-dominant singular values, those obtained in RM are larger than the HM ones. The corresponding time-domain right singular vectors also behave quite differently. For the RM case, the non-dominant right singular vectors are mainly due to clutter and not included in the approximation used for the excitation signal. Therefore, the first right singular vector can be used as the excitation signal with appropriate amplitudes and phase shifts dictated by the left singular vectors or TD-DORT eigenvectors at the desired frequency (Figure 8).

Figure 8. Phase distribution of the most significant TD-DORT eigenvector obtained at the central frequency (left) and those of the left singular vectors of the individual (middle) and full (right) SF-MDMs.

Figure 8. Phase distribution of the most significant TD-DORT eigenvector obtained at the central frequency (left) and those of the left singular vectors of the individual (middle) and full (right) SF-MDMs.

It is worth reminding that the excitation signals obtained via SF-MDM processing do not require phase smoothing, as opposed to the TD-DORT method. In certain cases where multiple scattering is very pronounced, phase smoothing might fail to yield coherent time-domain excitation signals, whereas SF-MDM processing will always yield coherent time domain signals. Figure 9 depicts the images obtained by TD-DORT and SF-imaging in homogeneous and random media. Note that these images correspond to the electric field distributions at the time of focusing and plotted in linear scale. Note that in terms of cross-range resolutions both methods provide similar performance. However, it was shown in [107] that in HM, TD-DORT slightly performs better. As alluded to above, in the RM case TD-DORT might suffer from phase-smoothing limitations under strong multiple scattering whereas SF-imaging is not affected by this limitation. We also note that it is not possible to claim that the SF-imaging provides improved resolution as the randomness (or multiple scattering) in the background increases. Especially, this is not possible as long as the approximate steering vector (Equation (14) is utilized. However, even for strong clutter cases where the TD-DORT might fail to generate coherent time-domain excitation signals, SF-imaging is able to generate coherent time-domain excitation signals which is an advantage of SF-imaging over the TD-DORT method. Finally, as long as the background medium is known and the exact background steering vectors are used for the imaging functionals, both methods might provide statistical stability thanks to the frequency decorrelation achieved by UWB operation.

It should be pointed out that some other signal space-based algorithms have also appeared in the literature, as in [122,123].

Figure 9. Images (in linear scale) obtained both in homogeneous (top row) and random (bottom row) media by TD-DORT (1st column), SF-imaging using individual (2nd column), and full SF-MDMs (3rd column).

Figure 9. Images (in linear scale) obtained both in homogeneous (top row) and random (bottom row) media by TD-DORT (1st column), SF-imaging using individual (2nd column), and full SF-MDMs (3rd column).

3.2. Null Space Methods

Both DORT and SF-imaging methods rely on the fact that, for well-resolved point-like scatterers, information on scatterer strength(s) and location(s) are partially encoded by the eigenvalues and associated eigenvectors in the TRO signal space [47]. Backpropagation of these eigenvectors yields the target images. However, the performance of these algorithms degrades if the well-resolvedness criterion is not met. In this case, the SS eigenvectors become linear combinations of the Green’s function vectors connecting scatterers to the TRA [40]. Backpropagation of such SS eigenvectors creates overlapped image fields which hamper target imaging and localization. On the other hand, regardless of the well-resolvedness criteria, the TRO NS is always orthogonal to the TRO SS, i.e., projection of any vector formed by the linear combination of SS eigenvectors onto the NS is (ideally) zero. This property is the basis of TR-based MUSIC methods, which provides better detection and localization properties than signal space based methods even for poorly-resolved scatterers (assuming homogeneous media). A mathematical summary of the TR-MUSIC is presented next.

3.2.1. TR-MUSIC

Recalling Equation (9), the TRO NS is formed by the eigenvectors having near zero eigenvalues as shown below:

$$NS\left(\omega \right)=\{{\mathbf{v}}_{M_s+1}\left(\omega \right),...,{\mathbf{v}}_N\left(\omega \right)\}\mathrm{with}{\lambda }_{M_s+1}^2\left(\omega \right)\approx ..\approx {\lambda }_N^2\left(\omega \right)\approx 0$$

(21)

Note that, an NS (and thus the possibility to utilize TR-MUSIC) exists as long as the number of significant eigenvalues is less than the number of TRA antennas, i.e., $M_s\le N$. This may not hold for non-isotropic scattering, where more than one eigenvalue may be associated with a single scatterer [49,85,86,124]. Similarly, with increasing clutter and/or noise, additional TRO eigenvalues may become significant, which makes it increasingly difficult to distinguish the clutter/noise contribution from that of the discrete scatterers. In this case, a threshold criteria can be set to determine the SS and NS. This threshold may depend on the frequency and clutter/noise level.

In order to form the scatterer images, TR-MUSIC method requires the knowledge of the exact background Green’s function vector $\mathbf{g}({\overline{X}}_s,\omega )$ (steering vector) at each search point ${\overline{X}}_s$ in the probed domain, which is defined by

$$\mathbf{g}({\overline{X}}_s,\omega )={\left[G({\overline{X}}_s,{\overline{R}}_1,\omega ),...,G({\overline{X}}_s,{\overline{R}}_N,\omega )\right]}^T$$

(22)

where $G(\overline{r},{\overline{r}}^{\prime },\omega )=|G(\overline{r},{\overline{r}}^{\prime },\omega )|e^{j\varphi (\overline{r},{\overline{r}}^{\prime },\omega )}$ is the Green’s function of the problem [16,57]. The conjugate of the steering vector provides the necessary phase and amplitude distribution for the array excitation to focus on the desired point. Since $\mathbf{g}({\overline{X}}_s,\omega )$ is often not known in an deterministic fashion in inhomogeneous media, approximate steering vectors ${\mathbf{g}}_0({\overline{X}}_s,\omega )$ can be employed instead.

As mentioned above, as the well-resolvedness criterion for the scatterers is weakened, SS eigenvectors become linear combinations of the Green’s function vectors connecting the scatterers to the TRA. This is illustrated in Figure 10.

Figure 10. Illustration of well-resolved and non-well-resolved cases for a scenario having $M_s=2$ scatterers and probed with $N=3$ antennas. While the SS is formed by the plane (P) formed by the first two eigenvectors (${\mathbf{v}}_1$ and ${\mathbf{v}}_2$), NS is formed by the eigenvector (${\mathbf{v}}_3$) orthogonal to P.

Figure 10. Illustration of well-resolved and non-well-resolved cases for a scenario having $M_s=2$ scatterers and probed with $N=3$ antennas. While the SS is formed by the plane (P) formed by the first two eigenvectors (${\mathbf{v}}_1$ and ${\mathbf{v}}_2$), NS is formed by the eigenvector (${\mathbf{v}}_3$) orthogonal to P.

Imaging using these eigenvectors creates an interference pattern that degrades scatterer(s) location(s) estimates. However, even for closely spaced scatterers, NS is still orthogonal to the SS. Therefore, inner products of the steering vectors with the NS eigenvectors would vanish only at the scatterer(s) location(s). This provides the TR-MUSIC imaging functional (pseudospectrum) as given below:

$$\mathbf{M}({\overline{x}}_s,\omega )={\left(\sum _{i=M_s+1}^N|〈{\mathbf{g}}_0({\overline{x}}_s,\omega ),{\mathbf{v}}_i^{*}\left(\omega \right)〉|\right)}^{-1}$$

(23)

Since this scheme uses a single (central) frequency (CF) ${\omega }_c$, the above is referred to as CF-TR-MUSIC, or simply CF-MUSIC in what follows.

In inhomogeneous media, ${\mathbf{g}}_0$ does not cancel out exactly the phase of the SS eigenvectors ${\mathbf{v}}_p(={\mathbf{g}}^{*}/\parallel \mathbf{g}\parallel$) nearby the scatterer(s) location(s) as $\mathbf{g}$ would. Similarly, the inner product of ${\mathbf{g}}_0$ with the NS eigenvectors may not necessarily produce a minimum (ideally zero) at the original scatterer locations. In addition, since $\mathbf{g}$ has phase factors that depend on the (random medium) realization, using ${\mathbf{g}}_0$ does not produce statistical stability. In order to achieve statistical stability, one needs to explore frequency decorrelation. This can be achieved by UWB operation and combining images obtained at different frequencies via the following frequency integration [40]:

$${\mathbf{M}}_{\Omega }\left({\overline{x}}_s\right)={\left({\int }_{\Omega }\sum _{i=M_s\left(\omega \right)+1}^N|〈{\mathbf{g}}_0({\overline{x}}_s,\omega ),{\mathbf{v}}_i^{*}\left(\omega \right)〉|d\omega \right)}^{-1}$$

(24)

where Ω is the bandwidth of operation. We call this strategy as UWB-TR-MUSIC, or simply UWB-MUSIC.

We apply both these TR-MUSIC algorithms to the same scenario of the previous subsection. Figure 11 shows the images obtained.

Figure 11. Images (in dB scale) obtained both in homogeneous (HM) and random media (RM) by CF-MUSIC and UWB-MUSIC.

Figure 11. Images (in dB scale) obtained both in homogeneous (HM) and random media (RM) by CF-MUSIC and UWB-MUSIC.

One of the first observations is that the dynamic range of the MUSIC images are larger than the signal space methods (Figure 9). This is due to the fact that MUSIC pseudospectrum utilizes the null subspace and the inner products around the scatterer locations should ideally be zero resulting in a very strong peak at the scatterer locations (especially for the homogeneous media). A more detailed comparison of the dynamic ranges for varying clutter conditions can be found in [108]. In a homogeneous medium case, both TR-MUSIC methods provide good co-range and cross-range resolutions and they overperform those of the TD-DORT and SF-imaging. In the random medium case, the clutter affects the image quality. CF-MUSIC fails to image the scatterer location due to multiple spurious peaks, which can be misinterpreted as discrete scatterer locations. Moreover, the images vary with different random medium realizations because the imaging functionals utilize ${\mathbf{g}}_0$ and not $\mathbf{g}$. Therefore, CF-MUSIC does not provide stable images under strong clutter. UWB-MUSIC, on the other hand, combines images obtained at different frequencies to yield a statistically stable image. This is produced at the expense of poorer co-range and cross-range resolutions (blurring) compared to the homogeneous case. Note that for both TR-MUSIC methods, increased clutter (i.e., decreased signal-to-clutter ratio) yields wider cross-range resolutions. Also, both TR-MUSIC methods fail to work in cases where NS is null, e.g., when all clutter eigenvalues are above the threshold. On the other hand, both DORT and SF-imaging methods provide similar co-range and cross-range information even for increasing clutter. As for processing times, it was shown in [108] that narrow-band methods (i.e., CF-MUSIC or CF-DORT) are faster than their UWB counterparts, and as long as the number SS eigenvectors is less than the number of NS eigenvectors, SS based methods are faster than the NS based methods. Additionally, TR-MUSIC is more stable against increasing dispersion and medium loss as compared to the SS based methods. Further details on the stability of the methods to perturbations can be found in [108,125]. There are also several other variant TR-MUSIC algorithms as considered, e.g., in [126,127]. Finally, in this paper, we have assumed that the number of distinct scatterers are less than the total number of antennas. For scenarios where the scatterers are more than the number of antennas, the methods explained in this paper may not be sufficient for imaging. As discussed in [72], a practical approach is to suppress the clutter first via an anti-focusing step and then attempt to image the scatterers of interest. This algorithm has shown to have good performance over conventional approaches when the number of scatterers is larger than the TRA sensors. So as a future study, we plan to incorporate this approach into the TR-imaging algorithms presented here. Especially, the effects of anti-focusing on the eigenspace and singular value/vectors structures will be focused.

4. Frequency Dispersion/Loss Compensation Techniques

The methods discussed so far have assumed that the intervening medium is lossless and stationary, hence the invariance of the wave equation under TR remains valid. However, dispersive and lossy media are often encountered in nature, such as in soils, rocks, ice, and most biological tissues. In such cases, the TR invariance is broken and conventional TR operation can not be directly applied. Several compensation methods (within certain limitations) have been discussed in the literature [22,75,128,129]. In this section, we review a more robust compensation method developed for UWB signals propagating in dispersive media [106,109].

A dispersive medium acts as a filter for UWB signals that propagate in it. For such a medium, the electric field at a single frequency is represented as

$$\overline{E}(\overline{r},\omega )=\overline{F}\left(\overline{r}\right)exp\left(-j\omega \sqrt{\mu ϵ\left(\omega \right)}\overline{r}\right)$$

(25)

where $\overline{F}\left(\overline{r}\right)$ is a frequency-independent amplitude coefficient and for simplicity the permittivity function $ϵ\left(\omega \right)$ is here assumed to be invariant over space ($ϵ(\overline{r},\omega )=ϵ\left(\omega \right)$) [130, pg. 18].

The non-zero imaginary part of $ϵ\left(\omega \right)$ yields additional attenuation on the received signals by the TRA relative to the non-dispersive case as shown in Figure 12. Moreover, the larger real part of $ϵ\left(\omega \right)$ causes an additional phase shift (delay). Since TR signals are phase conjugated coherently along the entire bandwidth, any additional phase shift induced from dispersion during the forward propagation is exactly compensated by the TR process. In a dispersive medium however, the signal is attenuated during both forward- and back-propagation. As a result, refocusing of the TR signals can be significantly degraded relative to the non-dispersive (lossless) case. To overcome this degradation, a compensation technique should be applied to act as an inverse filter with respect to the attenuation. Additionally, dispersive effects are also cumulative in time, i.e., the longer a signal travels in the medium, the more it is attenuated. In other words, a dispersive medium acts as a frequency-dependent, cumulative-in-time filter to the propagating signals. This requires the use of space and frequency dependent compensation filters. The implementation of the compensation method does not change the basic TR experiment, i.e., forward and backpropagation steps stay the same. But received and recorded signals after the forward propagation are properly modified (compensated) and then retransmitted. During the forward propagation of the transmitted short pulse in the dispersive medium, different multipaths go through different (dispersive) relative attenuations depending on the total distance (or time) traveled in the medium. In other words, received signals at later times (or equivalently, at further distances) need more strong “compensation” than those at earlier times. Ideally, different filters should be applied for each particular instant (sample), but this is impractical. Instead, time-windowing can be applied to received signals and different frequency-dependent filters used on each window. The combined time-windowing/filtering process can be best achieved after applying a short-time Fourier transform (STFT) [132] to the received signals at the TRA.

Figure 12. Magnitude of the exponential term $exp(-j\omega \sqrt{\mu ϵ\left(\omega \right)}\overline{r})$ plotted with respect to frequency and spatial distance for complex $ϵ\left(\omega \right)$ corresponding to the Puerto Rico soil with moisture level of %2.5 [131]. Space- and frequency-dependent attenuation is observed.

Figure 12. Magnitude of the exponential term $exp(-j\omega \sqrt{\mu ϵ\left(\omega \right)}\overline{r})$ plotted with respect to frequency and spatial distance for complex $ϵ\left(\omega \right)$ corresponding to the Puerto Rico soil with moisture level of %2.5 [131]. Space- and frequency-dependent attenuation is observed.

The whole compensation process is relatively simple and starts by taking the discrete STFT, $X[k,i]$, of the (sampled) received signal at each TRA antenna $x\left[n\right]$, which is given as

$$X[k,i]=\sum _{n}x\left[n\right]w_i\left[n\right]e^{-j(2\pi /L)nk}\mathrm{for}i=1,..,M$$

(26)

where $w_i\left[n\right]$ is the $i^{th}$ window, L is the FFT length, and M is the total number of windows. Note that M depends on the window length and overlapping factor (between windows). The product $x\left[n\right]w_i\left[n\right]=s_i\left[n\right]$ is referred to as the $i^{th}$ windowed signal. Each windowed signal has traveled a different distance amount to reach the receiver and each of its frequency components has gone through different attenuation. Therefore, each windowed signal has to be filtered with a different space and frequency dependent filter $H[k,i]$ obtained for the $i^{th}$ window. In the Fourier domain, the filtering is carried out by

$$Y[k,i]=X[k,i]H[k,i]\mathrm{for}i=1,..,M$$

(27)

The final, dispersion-compensated signal $x_c\left[n\right]$ is then obtained by an inverse STFT:

$$x_c\left[n\right]=\sum _{i=1}^M\frac{1}{w_i\left[n\right]}\sum _{k}Y[k,i]e^{j(2\pi /L)nk}$$

(28)

A block diagram of the whole procedure is given in Figure 13.

Filters play an important role here and to obtain them, the dispersive characteristics of the media should be (ideally) known or (practically) estimated. One way to obtain these filters is to compare the solution of the wave equation in a homogeneous reference medium having the exact dispersive model of interest with that in a nondispersive reference medium, as carried out in [109].

Figure 13. The block diagram of the dispersion compensation method.

Figure 13. The block diagram of the dispersion compensation method.

Now, let us illustrate the compensation procedure by considering the received signal at one of the TRA elements is considered in more detail. The received signal is windowed using five Hamming windows [132] of length $256\Delta t$ each and an overlapping factor of $0.5$, as shown in Figure 14.

Figure 14. (Left:) One of the original received signal and the employed Hamming windows with overlapping factor of 0.5; (right:) corresponding windowed signals. Note the amplitude difference for each windowed signal.

Figure 14. (Left:) One of the original received signal and the employed Hamming windows with overlapping factor of 0.5; (right:) corresponding windowed signals. Note the amplitude difference for each windowed signal.

For each windowed signal, a filter is designed as follows: First, for each central point of window $w_i\left[n\right]$, a corresponding effective distance is found for the specific dispersive characteristics. The ratio between the frequency domain representations of the signals received in the dispersive and non-dispersive reference media at this effective distance gives the frequency-dependent attenuation due to dispersion to be compensated by $H[k,i]$. Note that, to avoid noise contamination, the amplification factors can be smoothly set to unity for frequencies where the spectral density are below acceptable levels (high-frequency ends) depending on the particular application. Then, these filters are applied to corresponding windowed signals to yield the filtered windowed signals shown in Figure 15.

Figure 15. Time and Frequency domain representations for some of the windowed signals and their compensated counterparts after space and frequency dependent filtering.

Figure 15. Time and Frequency domain representations for some of the windowed signals and their compensated counterparts after space and frequency dependent filtering.

These filtered windowed signals are then inverse Short-Time Fourier transformed to obtain the final compensated signal, which is shown in Figure 16 along with the original signal received in the dispersive medium and a reference signal that would have been received in the non-dispersive case.

Figure 16. Signal received by one of the TRA antennas in dispersive medium, corresponding compensated signal and the reference signal that would be received in non-dispersive medium in (left) time and (right) frequency domains.

Figure 16. Signal received by one of the TRA antennas in dispersive medium, corresponding compensated signal and the reference signal that would be received in non-dispersive medium in (left) time and (right) frequency domains.

As seen in these plots, apart from the phase shift (which is automatically compensated by the TR process), the amplitude of the compensated signal is much closer to the reference signal amplitude than that of the original signal. Any discrepancy between the compensated and the reference signal is due to the finite window length used. A better agreement can be obtained with a larger number of windows, however, trade-off is the increased computational cost. This procedure is applied to all the signals received by the TRA. Note that, since backpropagation occurs in the same dispersive medium, additional attenuation is produced on the signals. Therefore, for each windowed signal, two different sets of compensation filters should be applied for each windowed signal. The first one compensates for propagation losses from the source (TRA) to scatterers and from scatterers back to the TRA; the second one compensates for propagation losses from TRA to scatterers. The compensation filters used for backpropagation should be different from those for the initial propagation, since the effective propagation duration for the backpropagation is approximately half of the initial propagation. The second set of filters for each windowed signal can be obtained by using half of the effective distances used for the first set of compensation filters. Once these compensated and TR signals are backpropagated into the same medium, the refocused fields around the scatterers have amplitudes much closer to those obtained in the nondispersive case, as shown in [106,109]. It should also be noted that compensation filters also amplify system noise; hence, they are of limited availability in scenarios where the medium dispersion reduces the signal to noise ratio below acceptable levels.

5. Conclusions

In this paper, we have provided a summary of some time-reversal techniques for UWB microwave remote sensing. For concreteness, we have focused most of the discussion on imaging scenarios consisting of obscured discrete targets in continuous random media, although TR techniques are applicable to many other remote sensing scenarios, as surveyed in the Introduction.

First, we considered the so-called “signal space” methods, viz., TD-DORT and SF-imaging, which respectively utilize space-space and space-frequency multistatic data matrices obtained from the TR array signals. For the TD-DORT method, the use of an eigenvector decomposition of the TR operator along the available bandwidth followed by a phase smoothing algorithm (on the eigenvectors associated with the significant eigenvalues) provides the required time-domain excitation signals to be used on the TRA for selective focusing of multiple targets. While TD-DORT performs well in relatively low clutter environments, it fails in environments with stronger clutter because of the limitations in the phase smoothing algorithm. SF-imaging, on the other hand, utilizes less conventional MDMs to provide coherent time-domain excitation signals regardless of the background clutter. Although the TD-DORT cross-range performance is slightly better than that provided by SF-imaging, the latter is favored in highly scattering environments.

Next, we considered the TR-MUSIC method as the prototypal “null space” based TR algorithm. In this case, the TR operator eigenvectors associated with near zero eigenvalues of the TR operator (as opposed to the significant eigenvalues, as in the “signal space” case) are utilized along with medium steering vectors to obtain desired target images. Under weak clutter conditions, TR-MUSIC outperforms TD-DORT and SF-imaging methods in terms of both co-range and cross-range resolutions. However, TR-MUSIC was found to be less stable under increased clutter.

Finally, we briefly discussed the use of compensation (inverse) filters as a means to extend the basic TR techniques to lossy media, where the TR invariance is broken.