Synthetic aperture imaging with intensity-only data

TL;DR

This paper introduces a novel algorithm for synthetic aperture imaging using intensity-only data, enabling coherent imaging by recovering phase information from intensity measurements across multiple frequencies and locations.

Contribution

It presents a new method to recover field cross-correlations and synchronize phase data from intensity-only measurements in synthetic aperture imaging.

Findings

Successfully recovers phase information up to a common phase

Enables coherent imaging with intensity-only data

Works across multiple frequencies and measurement locations

Abstract

We consider imaging the reflectivity of scatterers from intensity-only data recorded by a single moving transducer that both emits and receives signals, forming a synthetic aperture. By exploiting frequency illumination diversity, we obtain multiple intensity measurements at each location, from which we determine field cross-correlations using an appropriate phase controlled illumination strategy and the inner product polarization identity. The field cross-correlations obtained this way do not, however, provide all the missing phase information because they are determined up to a phase that depends on the receiver's location. The main result of this paper is an algorithm with which we recover the field cross-correlations up to a single phase that is common to all the data measured over the synthetic aperture, so all the data are synchronized. Thus, we can image coherently with data over all frequencies and measurement locations as if full phase information was recorded.

Figures

Equations

$${\bf{P}}^{j}=\left[P(\mbox{\boldmath{$x$}}_{j},\omega_{1})\quad P(\mbox{\boldmath{$x$}}_{j},\omega_{2})\quad\ldots\quad P(\mbox{\boldmath{$x$}}_{j},\omega_{S})\right]$$

$$\displaystyle[\mbox{\boldmath{$M$}}^{j}]_{ll^{\prime}}\equiv\displaystyle m^{j}_{ll^{\prime}}=\overline{P(\mbox{\boldmath{$x$}}_{j},\omega_{l})}P(\mbox{\boldmath{$x$}}_{j},\omega_{l^{\prime}}),$$

$$\displaystyle\mbox{Re}(m^{j}_{ll^{\prime}})=\frac{1}{2}\left(|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l+l^{\prime}}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l^{\prime}}|^{2}\right)\,,$$

$$\displaystyle\mbox{Im}(m^{j}_{ll^{\prime}})=\frac{1}{2}\left(|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l-{\imath}l^{\prime}}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l^{\prime}}|^{2}\right)\,.$$

$$\frac{m^{j}_{ll^{\prime}}}{|P(\mbox{\boldmath{$x$}}_{j},\omega_{l})|}=P(\mbox{\boldmath{$x$}}_{j},\omega_{l^{\prime}})\frac{\overline{P(\mbox{\boldmath{$x$}}_{j},\omega_{l})}}{|P(\mbox{\boldmath{$x$}}_{j},\omega_{l})|},\ l^{\prime}=1,\ldots,S.$$

$$b^{j}_{l}=P(\mbox{\boldmath{$x$}}_{j},\omega_{l})e^{\imath\theta_{j}},\mbox{ for }l=1,\ldots,S$$

$$P(\mbox{\boldmath{$x$}}_{j},\omega_{l})=\sum_{q=1}^{Q}\tilde{\rho}_{q}e^{\imath 2\frac{\omega_{l}}{c}r^{j}_{q}},$$

$$\widehat{{\rho }_j}({\kappa }_l)=\sum _{q=1}^Q{\tilde{\rho }}_q{e}^{ı{\kappa }_l{r}_q^j}.$$

$${\cal A}^{j}\tilde{{\mbox{\boldmath{$\rho$}}}}^{j}={\mbox{\boldmath{$b$}}}^{j},$$

$${\cal A}^{j}=\left[\begin{array}[]{cccc}e^{\imath 2\frac{\omega_{1}}{c}r^{j}_{1}}&e^{\imath 2\frac{\omega_{1}}{c}r^{j}_{2}}&\cdots&e^{\imath 2\frac{\omega_{1}}{c}r^{j}_{Q}}\\ e^{\imath 2\frac{\omega_{2}}{c}r^{j}_{1}}&e^{\imath 2\frac{\omega_{2}}{c}r^{j}_{2}}&\cdots&e^{\imath 2\frac{\omega_{2}}{c}r^{j}_{Q}}\\ \vdots&&&\vdots\\ e^{\imath 2\frac{\omega_{S}}{c}r^{j}_{1}}&e^{\imath 2\frac{\omega_{S}}{c}r^{j}_{2}}&\cdots&e^{\imath 2\frac{\omega_{S}}{c}r^{j}_{Q}}\\ \end{array}\right],$$

$${R}_{k}^{j}=|\mbox{\boldmath{$y$}}_{k}-\mbox{\boldmath{$x$}}_{j}|.$$

$$\begin{array}[]{l}\displaystyle\mbox{set }i=1\mbox{ and }r_{i}^{j}=R_{i}^{j}\\ \mbox{for k=2 to K}\\ \hskip 14.22636pt\mbox{if }|R_{k}^{j}-R_{k-1}^{j}|>\epsilon\\ \hskip 19.91684pt\mbox{set }i=i+1\mbox{ and }r_{i}^{j}=R_{k}^{j}\\ \hskip 14.22636pt\mbox{end}\\ \mbox{end}\end{array}$$

$$\sum_{q=1}^{Q}\tilde{\rho}^{j}_{q}\approx e^{\imath\theta_{j}}\frac{1}{h_{c}}\int_{IW}{\mbox{\boldmath{$\rho$}}}_{0}d\vec{\bf y},$$

$$c_j=\frac{{\sum }_{q=1}^Q{\tilde{\rho }}_q^j}{{\sum }_{q=1}^Q{\tilde{\rho }}_q^1}\stackrel{\text{eqref}eq:Arhoj}{=}{e}^{ı({\theta }_j-{\theta }_1)},j=1,\dots ,N,$$

$$\overline{c_{j}}b^{j}_{l}=P(\mbox{\boldmath{$x$}}_{j},\omega_{l})e^{\imath\theta_{1}},\,\,\forall\,j=2,\ldots,N\,\mbox{and}\,\,\forall\,\,l=1,\ldots,S.$$

$$\begin{array}[]{ll}\displaystyle P^{h}(\mbox{\boldmath{$x$}}_{1},\omega_{l})&\displaystyle=b^{1}_{l},\,\,\forall\,\,l=1,\ldots,S.\\ \displaystyle P^{h}(\mbox{\boldmath{$x$}}_{j},\omega_{l})&\displaystyle=\overline{c_{j}}b^{j}_{l},\,\,\forall\,j=2,\ldots,N\,\mbox{and}\,\,l=1,\ldots,S.\end{array}$$

$${\rho}^{KM}({\mbox{\boldmath{$y$}}}_{k})=\sum_{j=1}^{N}\sum_{l=1}^{S}e^{-\imath 2\frac{\omega_{l}}{c}|\mbox{\boldmath{$x$}}_{j}-\mbox{\boldmath{$y$}}_{k}|}P^{h}(\mbox{\boldmath{$x$}}_{j},\omega_{l}),$$

$${\cal A}\mbox{\boldmath$\rho$}=\mbox{\boldmath$b$}\,,$$

$${\cal A}=\left[\begin{array}[]{cccc}e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{1}|}&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{2}|}&\cdots&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{K}|}\\ e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{2}-\mbox{\boldmath{$y$}}_{1}|}&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{2}-\mbox{\boldmath{$y$}}_{2}|}&\cdots&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{2}-\mbox{\boldmath{$y$}}_{K}|}\\ \vdots&\vdots&&\vdots\\ e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{1}|}&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{2}|}&\cdots&e^{\imath 2\frac{\omega_{1}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{K}|}\\ e^{\imath 2\frac{\omega_{2}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{1}|}&e^{\imath 2\frac{\omega_{2}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{2}|}&\cdots&e^{\imath 2\frac{\omega_{2}}{c}|\mbox{\boldmath{$x$}}_{1}-\mbox{\boldmath{$y$}}_{K}|}\\ \vdots&\vdots&&\vdots\\ e^{\imath 2\frac{\omega_{S}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{1}|}&e^{\imath 2\frac{\omega_{S}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{2}|}&\cdots&e^{\imath 2\frac{\omega_{S}}{c}|\mbox{\boldmath{$x$}}_{N}-\mbox{\boldmath{$y$}}_{K}|}\\ \end{array}\right]$$

$$b_{i}=b_{(l-1)N+j}=P^{h}(\mbox{\boldmath{$x$}}_{j},\omega_{l}),\,\,j=1,\ldots,N,l=1,\ldots,S.$$

$$\min\|\mbox{\boldmath$\rho$}\|_{\ell_{1}}\quad\text{subject to }\quad{\cal A}\mbox{\boldmath$\rho$}=\mbox{\boldmath{$b$}}\,,$$

$${\eta }_{\tau }(y_i)=\text{sign}(y_i)\max \{0,|y_i|-\tau \}.$$

$$[\mbox{\boldmath{$M$}}^{j}]_{ll^{\prime}}\equiv\displaystyle m^{j}_{ll^{\prime}}=\overline{P(\mbox{\boldmath{$x$}}_{j},\omega_{l})}P(\mbox{\boldmath{$x$}}_{j},\omega_{l^{\prime}}),\quad l,l^{\prime}=1,\dots,S\,\,,$$

$$\mbox{Re}(m^{j}_{ll^{\prime}})=\frac{1}{2}\left(|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l+l^{\prime}}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l^{\prime}}|^{2}\right)$$

$$\mbox{Im}(m^{j}_{ll^{\prime}})=\frac{1}{2}\left(|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l-{\imath}l^{\prime}}|^{2}-|{\bf{P}}^{j}\cdot\mbox{\boldmath{$e$}}_{l}|^{2}-|{\bf{P}}^{j}\mbox{\boldmath{$e$}}_{l^{\prime}}|^{2}\right).$$

$$m_{l{l}^{\prime }}^j=\frac{m_{l1}^j{m}_{1l^{\prime }}^j}{m_{11}^j},$$

$$I_0$$

$$I_{\pi /2}$$

$$\left[\begin{array}[]{ll}\alpha_{1}&-\alpha_{2}\\ \alpha_{2}&\alpha_{1}\end{array}\right]\left[\begin{array}[]{l}\cos(\Delta\phi^{j}_{ll^{\prime}})\\ \sin(\Delta\phi^{j}_{ll^{\prime}})\end{array}\right]=\left[\begin{array}[]{l}\beta_{1}\\ \beta_{2}\end{array}\right]$$

$${\alpha }_1={\int }_{t_0}^{t_0+\Delta t}\cos (({\omega }_l-{\omega }_{l^{\prime }})t)dt,{\alpha }_2={\int }_{t_0}^{t_0+\Delta t}\sin (({\omega }_l-{\omega }_{l^{\prime }})t)dt$$

$$\beta_{1}=\int_{t_{0}}^{t_{0}+\Delta t}\frac{I_{0}-\left|P(\mbox{\boldmath{$x$}}_{j},\omega_{l})\right|^{2}-\left|P(\mbox{\boldmath{$x$}}_{j},\omega_{l^{\prime}})\right|^{2}}{2\left|P(\mbox{\boldmath{$x$}}_{j},\omega_{l})\right|\left|P(\mbox{\boldmath{$x$}}_{j},\omega_{l^{\prime}})\right|}dt$$

Full text

1 Introduction

We consider a multifrequency phaseless synthetic aperture imaging system composed of a single transmitter/receiver element that operates at microwave frequencies. The system only records the intensities of the signals, and forms the images by combining the data coherently over the entire synthetic aperture. With the proposed computational imaging approach, we show that using intensity-only data is as good as coherent imaging with full (phase and amplitude) data; cross-range and range resolutions are still inversely proportional to the synthetic aperture size and to the bandwidth of the recorded signals, respectively.

Imaging with microwaves is of interest in security applications because they can penetrate through materials that are opaque at visible wavelengths, allowing the detection of concealed objects under cloths and inside luggage [21, 27]. Microwaves have also been successfully used for through-wall imaging [30], breast cancer detection in medical imaging [25, 16], and space surveillance [8]. An arbitrary number of transmitters and receivers is used in general, typically arranged in a planar or spherical geometry, while a moving single transmitter/receiver element can also be used so as to form a synthetic aperture. This last mode records less information but it can still provide data with enough frequency and space diversity for imaging.

In all applications, control and detection of phases is essential to image coherently for good resolution. However, to maintain phase coherence during the whole data acquisition process may be difficult or impossible. This is the case, for example, when high frequency signals, above $30$ GHz, and large synthetic apertures are used to form high quality images [18, 32]. Other situations in which phase measurements may not be reliable arise when there is uncertainty in the antenna location or in the signal trajectory.

The conventional way to overcome lack of coherent phase measurements is to resort to computational imaging systems that use intensity-only data recorded with simpler system architectures. The missing phases are then reconstructed with the well-known phase retrieval algorithms, which rely in an essential way on prior information about the object to be imaged. This is the case in various fields including crystallography, optical imaging, astronomy, and electron microscopy.

We recount briefly the basic facts about phase retrieval. In this problem, the objective is to reconstruct a signal from its power spectrum. This is, of course, an ill-posed problem because there is not enough information to recover the true signal. To resolve this issue, one can invoke prior knowledge about the signal and look for a solution by using a phase retrieval algorithm. The most widely used algorithms are the alternating projection algorithms introduced by Gerchberg and Saxton [15] and by Fienup [11, 12] that project the iterates on intensity data (or prior information) sequentially in both the real and the Fourier spaces. Although these algorithms are efficient and flexible for reconstructing the missing phases in the data, and performance is often good in practice, they may not converge or even get close to the true, missing phases because the problem is non-convex.

In [29], the authors use the quadratic approach which formulates the phase retrieval problem as a non-linear (quadratic) inverse problem with data the square amplitude of the near field. For this quadratic inverse problem, the presence of local minima can be avoided by increasing the number of independent data (cf. [29] and references therein). A convex, non-iterative approach that guarantees exact recovery in the case of sparse reflectivity was proposed in [6, 4], but its computational cost is high when the problem is large.

Holographic methods can also be used. These are interferometric approaches that obtain phase field differences from intensity data. For example, phase information can be recovered by superposing the scattered signals with a known and well-controlled reference signal. Holography was first introduced to increase the resolution of electron microscopes [14]. It was a two-step method for recording the phase of optical signals, since photographic film is not sensitive to complex amplitudes but to intensities. Later, holography was used with microwaves for measuring antenna radiation patterns [19, 1]. For recent adaptations of microwave holographic techniques, we refer to [7, 17, 28, 18].

On the other hand, the Wirtinger Flow (WF) phase retrieval algorithm of [5] is used in [13] for microwave imaging with a frequency-diverse metasurface antenna. The antenna produces spatially diverse radiation patterns that vary as a function of the frequency sampled over the operational K-band (17.5-26.5 GHz). In [32], the authors use the more recent sparse WF algorithm proposed in [31] that allows to reduce the computational cost of the method.

We follow here a different strategy inspired by interferometry. The key idea is to resolve the non-uniqueness of the phase retrieval problem creating redundancy in the data by illuminating the image multiple times. Indeed, by using an appropriate illumination strategy and the inner product polarization identity, the missing phase information can be uniquely determined, up to a global phase [26, 23, 24]. The polarization identity is a well known formula in mathematics that relates the inner product of two vectors with their norms. It is not related to the polarization of electromagnetic waves. Special forms of the polarization identity have also been used in [20] where antenna phase patterns are obtained from the responses to two probe antennas recorded at three power detectors, and in [7] where intensity data are collected with two probes at fixed offset moving over an arbitrary scanning surface.

In this paper, we present a new computational imaging approach to accurately reconstruct the reflectivity of scatterers with synthetic aperture, intensity-only data. The method has two stages. First, from the intensity data at each source-detector position, we recover field cross-correlations corresponding to coherent sources of different frequencies. This is achieved using a special sequence of illuminations that exploits the frequency diversity available on the transmitter side. The recovered field cross-correlations are the same as the ones obtained from full data, up to a phase that is different at each receiver location. Hence, at this point, these field cross-correlations cannot be combined coherently to image the reflectivity. To use them coherently over all the synthetic aperture they need to be synchronized or aligned first.

This is the second step of the proposed method in which all the phases that depend on the receiver locations are referred to a single global phase. The main idea of the second step is to refer the total reflectivity estimated at each receiver location to the total true reflectivity, which is a common quantity for all measurement locations. This is the main contribution of the paper. With the strategy proposed here, we show that imaging with intensity-only data is as good as imaging with full-phase data.

The paper is organized as follows. In Section 2, we explain the proposed method to obtain coherent cross-correlations when one element that transmits multifrequency microwaves and measures only intensities is used to collect the data on a synthetic aperture. In Section 3, we discuss two imaging methods traditionally used when full (phase and amplitude) data are available. We use these methods with the recovered cross-correlations. In Section 4, we show the results of our numerical experiments. Section 5 contains our conclusions.

2 Multi-frequency interferometric synthetic aperture imaging

In the next two subsections we present in detail the two stages of the proposed approach. The goal is to determine the reflectivity $\rho$ within a region of interest, called the image window IW, from multiple intensity-only measurements at different locations \mbox{\boldmath{x}}_{j}, $j=1,\ldots,N$, and frequencies $\omega_{l}$, $l=1,\ldots,S$, with a total number of $D=N\cdot S$ data (see Figure 1). We denote by |P(\mbox{\boldmath{x}}_{j},\omega_{l})| the amplitude of the signal received at location \mbox{\boldmath{x}}_{j} when a signal of frequency $\omega_{l}$, unit amplitude, and zero phase is emitted from the same location.

For imaging purposes, the IW is discretized using a uniform grid of $K$ points {\mbox{\boldmath{y}}}_{k}, $k=1,\ldots,K$. The unknown is the reflectivity vector \mbox{\boldmath{\rho}}=[\rho_{1},\ldots,\rho_{K}]^{t}\in\mathbb{C}^{K}, whose entries are the values of the reflectivity \rho_{k}=\rho({\mbox{\boldmath{y}}}_{k}) on the grid points {\mbox{\boldmath{y}}}_{k}, $k=1,\ldots,K$. We assume that $K>D$, and often we have $K\gg D$. Moreover, we assume that the unknown reflectivity vector is $M$-sparse with $M\ll K$. This is often true in applications where the reflectivity to be imaged does not occupy the entire scene but rather a small part of the IW.

2.1 Multi-frequency field cross-correlations

We pursue here the idea of [26, 23, 24], where it is shown that field cross-correlations can be obtained from intensity-only measurements by using an appropriate protocol of illuminations and the polarization identity. In these works, the illumination strategy was implemented with an array in which all its elements were used to emit and receive signals.

In synthetic aperture imaging, however, there is an inherent loss of possible illuminations because only one transmitter/receiver element is used. This lack of flexibility, or diversity, in illuminations might at first suggest that the data cannot be used coherently to form images when only intensities are recorded.

To consider this issue further, we introduce the row vector

[TABLE]

with $S$ components. The entry P(\mbox{\boldmath{x}}_{j},\omega_{l}) corresponds to the signal recorded at \mbox{\boldmath{x}}_{j}, including phases, when a unit amplitude and zero phase signal of frequency $\omega_{l}$ is sent from the same location \mbox{\boldmath{x}}_{j}. The full phase retrieval problem consists on determining, at once, the $N\times S$ phases of the phaseless measurements \left|P(\mbox{\boldmath{x}}_{j},\omega_{l})\right| for the $N$ measurement positions and the $S$ frequencies.

In this paper, we use a different strategy. First, we recover the field cross-correlation matrices at each receiver, up to a phase that depends on the receiver location but not on the frequency. In a second step, all these phases are referred to a single one; and this allows for coherent imaging.

Specifically, the first step of the proposed methodology consists on recovering, at every fixed location \mbox{\boldmath{x}}_{j}, the field cross-correlation matrix

[TABLE]

for $l,l^{\prime}=1,\dots,S$. To do so, we use the diversity over illuminations with different frequencies, and the inner product polarization identity

[TABLE]

[TABLE]

Here, \mbox{\boldmath{e}}_{l}\in\mathbb{C}^{S} is the vector with $1$ in the $l$-th coordinate and [math]’s elsewhere. It represents a signal of unit amplitude and zero phase at frequency $\omega_{l}$. In (3)-(4), ${\imath}=\sqrt{-1}$, \mbox{\boldmath{e}}_{l+l^{\prime}}=\mbox{\boldmath{e}}_{l}+\mbox{\boldmath{e}}_{l^{\prime}}, and \mbox{\boldmath{e}}_{l-{\imath}l^{\prime}}=\mbox{\boldmath{e}}_{l}-{\imath}\mbox{\boldmath{e}}_{l^{\prime}}. The vector \mbox{\boldmath{e}}_{l+l^{\prime}}=\mbox{\boldmath{e}}_{l}+\mbox{\boldmath{e}}_{l^{\prime}} refers to sending simultaneously signals of unit amplitude and zero phase at two frequencies $\omega_{l}$ and $\omega_{l^{\prime}}$, while the vector -{\imath}\mbox{\boldmath{e}}_{l^{\prime}} denotes a $-\pi/2$ phase shift in the signal of frequency $\omega_{l^{\prime}}$.

When two signals of distinct frequencies $\omega_{l}$ and $\omega_{l^{\prime}}$ are sent simultaneously to probe the medium, the intensity at the receivers oscillates at the difference frequency $\omega_{l}-\omega_{l^{\prime}}$. Therefore, our imaging system should be able to resolve intensities that oscillate at frequencies of the order of the total bandwidth $B$, which determines the range resolution $c/B$ of the imaging system, where $c$ is the signal speed. In Appendix B, we describe in detail the measurement process with which we recover all the elements of the cross-correlation matrix \mbox{\boldmath{M}}^{j}. We do not assume, however, cross-correlation measurements between different positions \mbox{\boldmath{x}}_{j} and, hence, we recover \mbox{\boldmath{M}}^{j} up to an unknown global phase $\theta_{j}$, i.e., a factor of the form $e^{i\theta_{j}}$, with $\theta_{j}$ depending on \mbox{\boldmath{x}}_{j}.

This means that, in principle, we can only recover asynchronized signals between the measurement locations. We remark that synchronization is needed regardless of whether the elements in \mbox{\boldmath{M}}^{j} are obtained through (3)-(4). It is a problem that arises due to the lack of information between signals sent and received at different positions. Hence, the cross-correlations over frequencies need to be synchronized or aligned over the receiver positions to image coherently.

This is the main difficulty in synthetic aperture imaging when the phases are not measured. We explain in the next subsection how to synchronize the signals in the frequency domain, i.e., how to recover $N$ individual phases from phase difference measurements of the form (2). Once all the signals are synchronized, the imaging problem is trivial as we show in Section 3.

We describe in Appendix A an illumination strategy that requires $3S-2$ illuminations to recover all the entries in \mbox{\boldmath{M}}^{j}. This is the minimum number of measurements needed per receiver location.

2.2 Location-dependent phase recovery

As discussed above, we can recover, up to a global phase that depends on the receiver location \mbox{\boldmath{x}}_{j}, field cross-correlations of the form (2) using the frequency diversity in the illuminations. Since the amplitudes |P(\mbox{\boldmath{x}}_{j},\omega_{l})| are known at every location \mbox{\boldmath{x}}_{j} for all the frequencies $\omega_{l}$, we can compute

[TABLE]

This means that full data can be recovered at each measurement location \mbox{\boldmath{x}}_{j} up to a global phase $\theta_{j}$ which is unknown. In other words,

[TABLE]

is known to us. To image coherently, we have to refer all the unknown phases $\theta_{j}$ to a single location. This is a synchronization problem with $N$ unknown phases instead of $N\times S$ unknown phases as for the full phase retrieval problem. We next explain how to solve this problem.

Remark 1

We assume here that the measurements |P(\mbox{\boldmath{x}}_{j},\omega_{l})| are kept only when they are above some threshold determined by the noise level. Otherwise, they are not used. In the numerical simulations discussed below it was not necessary to discard any measurements, even in the [math]dB case.

By linearizing the scattering problem and assuming that multiple scattering is negligible, the data is given by

[TABLE]

where $c$ is the velocity in a homogeneous medium, and $\tilde{\rho}_{q}$ is the integral of the reflectivity on the sphere of radius $r^{j}_{q}$ centered at \mbox{\boldmath{x}}_{j}; see Figure 2. Typically for sparse $\rho$, $\tilde{\rho}_{q}$ is a reflectivity of a single scatterer. We assume here that the reflectivity is frequency independent.

It follows from (6) that the data P(\mbox{\boldmath{x}}_{j},\omega_{l}) is the Fourier coefficient of the reflectivity \tilde{{\mbox{\boldmath{\rho}}}} corresponding to wavenumber $\kappa_{l}=2\frac{\omega_{l}}{c}$, i.e.,

[TABLE]

Therefore, at each source-receiver position we have to solve a one dimensional problem to recover $\tilde{\rho}_{q}$, $q=1,\ldots,Q$, from the processed data (5). Again, these are data with phases that are well defined up to global phases $\theta_{j}$ that are not known. This means that the data are trains of asynchronous spikes, each train corresponding to a measurement location \mbox{\boldmath{x}}_{j}.

Still, we can determine the vector \tilde{\mbox{\boldmath\rho}}=[\tilde{\rho}_{1},\tilde{\rho}_{2},\dots,\tilde{\rho}_{Q}] by solving the following linear system

[TABLE]

with sensing matrix

[TABLE]

and data vector {\mbox{\boldmath{b}}}^{j} with components given by (5). The superscript $j$ is used to stress that the linear systems (7) uses data recovered at location \mbox{\boldmath{x}}_{j}.

These linear systems are underdetermined and, hence, there are infinitely many scatterer’s configurations that match the data. However, if the true reflectivity \mbox{\boldmath\rho}_{0} is sparse, with only a few components different than zero, an $\ell_{1}$-minimization approach can find their unique sparse solution. Exact recovery is guaranteed under the assumption that the mutual coherence111The mutual coherence of ${\cal A}$ is defined as \max_{i\neq j}|\langle\mbox{\boldmath{a}}_{i},\mbox{\boldmath{a}}_{j}\rangle| with \mbox{\boldmath{a}}_{i}\in\mathbb{C}^{N} the columns of ${\cal A}$ normalized to one, so that \|\mbox{\boldmath{a}}_{i}\|_{\ell_{2}}=1 $\forall\,i=1,\ldots,K$. of the matrices ${\cal A}^{j}$ are smaller than $1/(2M)$, where $M$ is the number of non zero components in a vector \tilde{{\mbox{\boldmath{\rho}}}}^{j}. For more details about $\ell_{1}$-minimization methods we refer to [9, 3, 10, 2]. In the simulations shown below, we use a generalized lagrangian multiplier algorithm (GeLMA) [22], described in Algortihm 1, to find the sparsest solution to (7).

The matrix ${\cal A}^{j}$ defined in (8) depends on the radii $r^{j}_{i}$, $i=1,\ldots,Q$, which are computed in the following way. Given an IW with discretization points \mbox{\boldmath{y}}_{k}, $k=1,\ldots,K$, we compute the distances from all points \mbox{\boldmath{y}}_{k} to the receiver location \mbox{\boldmath{x}}_{j},

[TABLE]

These form the components of a vector in $\mathbb{R}^{K}$. We then sort the components of this vector in ascending order and keep only the entries that appear with multiplicity larger than one. In practice, we only keep the entries that differ from each other by at least a level $\epsilon$. The value of $\epsilon$ should be small enough so we do not disregard many components since that would affect the accuracy of the reconstruction, but it cannot be very close to zero because we do not want the columns of ${\cal A}^{j}$ to be almost parallel. Note that $\epsilon$ has units of length and should be chosen to be small with respect to the wavelength and the pixel size, so that neglecting distances that differ less than $\epsilon$ does not affect the accuracy of the recovered phases. Assuming the $R_{k}^{j}$ are ordered then this can be done as follows,

[TABLE]

This process generates the radii $r_{i}^{j}$, $i=1,\ldots,Q$, with $Q\leq K$, that depend on the receiver locations \mbox{\boldmath{x}}_{j}, $j=1,\ldots,N$.

Once the solution vectors \tilde{{\mbox{\boldmath{\rho}}}}^{j} have been found, we compute the total reflectivity within the IW by summing all the components of the vectors \tilde{{\mbox{\boldmath{\rho}}}}^{j}. That is, for each receiver location we compute the scalar

[TABLE]

with $h_{c}$ a constant that depends on the discretization. The key point here is that for all receiver positions we can compute an approximation to the total reflectivity \int_{IW}{\mbox{\boldmath{\rho}}}_{0}d\vec{\bf y}, up to unknown phase factors $e^{\imath\theta_{j}}$, $j=1,\ldots,N$.

Thus, we can refer all the recovered quantities (5) to a same global phase with no physical meaning for imaging purposes. Indeed, let us define the quantities

[TABLE]

by dividing the total reflectivities associated to every location \mbox{\boldmath{x}}_{j} by the total reflectivity obtained from the measurements recorded at \mbox{\boldmath{x}}_{1}. The choice of $j=1$ in the denominator in (12) is, of course, arbitrary. With this choice, $c_{1}=1$. Then, by multiplying the recovered data (5) by $\overline{c_{j}}$ we get

[TABLE]

This defines the holographic data

[TABLE]

The phases in (14) are now coherent over different receiver positions and frequencies! Thus, the unknown reflectivity $\rho$ can be reconstructed as if data with phases were recorded.

We want to emphasize that the proposed methodology allows one to produce holographic data from intensity measurements. This is of considerable importance since: i) intensity data are much easier to obtain, and can be recorded with less expensive equipment (sensors) than data obtained with holographic techniques, ii) holographic data contain coherent phase information and allow us to obtain depth resolved reconstructions, and iii) the proposed methodology does not need any prior information about the sought reflectivity. We compare next the performance of different imaging methods using (14) as data.

3 Full phase synthetic aperture imaging methods

Once the holographic data (14) are obtained, the unknown reflectivity can be reconstructed with any imaging method as if the data with phases were recorded with a synthetic array aperture. Here we show results obtained with the frequently used Kirchhoff migration (KM) imaging method and the $\ell_{1}$-optimization approach.

KM is a direct imaging method which can be written as

[TABLE]

where |\mbox{\boldmath{x}}_{j}-\mbox{\boldmath{y}}_{k}| is the distance between the measurement location \mbox{\boldmath{x}}_{j} and the search point \mbox{\boldmath{y}}_{k} in the IW. The image {\mbox{\boldmath\rho}}^{KM}=[{\rho}^{KM}({\mbox{\boldmath{y}}}_{1}),{\rho}^{KM}({\mbox{\boldmath{y}}}_{2})\dots,{\rho}^{KM}({\mbox{\boldmath{y}}}_{K})] is an approximation to the true reflectivity vector \mbox{\boldmath\rho}_{0}.

We also form an image by promoting a sparse solution to the linear system

[TABLE]

where \mbox{\boldmath\rho}\in\mathbb{C}^{K} is the sought reflectivity vector,

[TABLE]

is the model matrix, and \mbox{\boldmathb}\in\mathbb{C}^{N\cdot S} is the recovered data vector whose components are

[TABLE]

We note that the KM solution (15) can also be written as \mbox{\boldmath\rho}^{KM}={\cal A}^{*}\mbox{\boldmathb}, where ${\cal A}^{*}$ is the complex conjugate transpose of ${\cal A}$.

To find the sparsest solution to the system (16), we solve the $\ell_{1}$-minimization problem

[TABLE]

using GeLMA described in Algorithm 1. This algorithm involves matrix-vector multiplications followed by a shrinkage-thresholding step defined by the operator

[TABLE]

GeLMA converges to the solution of (19) independently of the regularization parameter $\tau$, see [22].

As discussed previously, the solution of (19) agrees with the exact $M$-sparse solution $\rho$ for noiseless data if the mutual coherence of the matrix ${\cal A}$ is smaller than $1/(2M)$.

4 Numerical simulations

We consider a high frequency microwave scanning regime with central frequency $f_{0}=50$GHz which corresponds to $\lambda_{0}=6$mm. We make measurements for $S=41$ equispaced frequencies covering a total bandwidth of $10$GHz using a single transmitter/receiver that is moving along a linear trajectory. The synthetic aperture is $a=20$cm, and the distance from its center to the center of the IW is $L=1$m; see Figure 3. We assume that the medium between the synthetic array and the IW is homogeneous. The size of the IW is $48\mbox{cm}\times 48\mbox{cm}$, and the pixel size is $6\mbox{mm}\times 6\mbox{mm}$. The measurements are gathered at $N=41$ equispaced locations. These parameters are typical in microwave scanning technology [18].

The numerical simulations are done in 3D, but with a 2D reflector geometry and a SAR trajectory on the plane of the scatterers as illustrated in Figure 3. In this configuration, the horizontal axis shows the range or depth of the scatterers, and the vertical axis shows their cross-range. Because the phases of the signals carry the information from the scatterers’ range, it is usually assumed that range cannot be determined from phaseless data. We see, however, that with the proposed computational imaging method both range and cross-range are obtained as if full phase information was available.

We assume that time-resolved intensities can be measured for a set of pulses of the form $\exp(\imath\omega_{l}t)\exp(-t^{2}/(2\sigma_{t})^{2})$, $\omega_{l}=\omega_{1},\omega_{2},\dots,\omega_{S}$, where $\sigma_{t}$ is the pulsewidth that is inversely proportional to the available bandwidth $B$. As discussed in Appendix A, to retrieve the phases for the $S$ frequencies we need to measure the intensities for $3S-2$ illuminations. A simple strategy consists of using a sufficient delay between successive illuminations so that the corresponding echoes are non-overlapping. Given an illumination at time [math], an estimate for the starting time for the scattered signal is $\frac{2L}{c}$, while its duration is of the order $\frac{2\,{\rm IW}_{\rm size}}{c}$, where ${\rm IW}_{\rm size}$ is the size of the IW. Thus, an estimate of the acquisition time per location is $(3S-2)\left(\frac{1}{2B}+\frac{2\,{\rm IW}_{\rm size}}{c}+\frac{2L}{c}\right)$. For the specific parameters used in the simulations this is about $(3S-2)10{\rm ns}=1.2\mu{\rm s}$.

In our simulations, we add to the signals that arrive to the receiver, including phases, mean zero gaussian noise corresponding to a SNR $=10$dB. Then, their intensities are computed and these are the data from which we form the images

Indeed, following the methodology described in Section 2, we recover from these intensities the holographic data (14), which have phases that are coherent over frequency and measurement locations. In our numerical examples, we used $\epsilon=0.001\mu{\rm m}$ in (10). Note that $\epsilon$ is small with respect to the wavelength and the pixel size, so that neglecting distances that differ less than $\epsilon$ does not affect the accuracy of the recovered phases. The Algortihm 1 is used with \tau=20\left<\left|{\cal A}^{*}\mbox{\boldmath{b}}\right|\right> and $\beta=\frac{1}{2\|{\cal A}\|_{2}}$. Here, ${\cal A}^{*}$ is the complex conjugate transpose of ${\cal A}$ and $<\cdot>$ denotes the mean. The termination criterium is \frac{\|\mbox{\boldmath{z}}_{k}-\mbox{\boldmath{z}}_{k-1}\|_{2}}{\|\mbox{\boldmath{z}}_{k-1}\|_{2}}\leq 1.0e-13, with \mbox{\boldmath{z}}_{k} and \mbox{\boldmath{z}}_{k-1} denoting the vector $z$ during the current and the previous iteration.

The results are shown in Figure 4. The top row of Figure 4 shows the distribution of targets we seek to find. The bottom left panel is the KM image, and the right panel is the image obtained with the $\ell_{1}$-minimization algorithm. As expected, KM shows resolution $\lambda_{0}L/a$ in the cross-range direction and $c/B$ in the range direction, which for our imaging setup corresponds to a resolution of $5\lambda_{0}$ in both directions. On the other hand, the image obtained with the $\ell_{1}$-minimization algorithm is almost exact.

To show the robustness of the proposed method when phases are not recorded, we consider next the same imaging configuration but assuming that the scatterers are displaced with respect to the grid points of the IW. Note that misplacements with respect to the grid amounts to a systematic modelling error that affects the accuracy of the recovered phases and, thus, deteriorates the image reconstruction. This is due to errors in the computed distances to the grid points ${R}_{k}^{j}$ (see (9)) used in the definition of the radii $r^{j}_{i}$ in the model matrices (8) and (17).

In Figure 5, we show results for scatterers that are displaced by $\lambda/8=0.75$mm (top row) and by $\lambda/2=3$mm (bottom row) with respect to the grid points in both range and cross-range directions. As expected, the reconstructions deteriorate as the displacement with respect to the grid increases, but they remain quite accurate even for the largest possible displacement value of half the grid size; see the bottom row plots in Figure 5.

Finally, we also add to the signals that arrive to the receiver corresponding to the displacement $\lambda/2=3$mm, including phases, mean zero gaussian noise corresponding to a SNR $=0$dB. The results shown in Figure 6 are also very good. They illustrate the robustness of the proposed methodology with respect to both additive noise and off-grid displacement errors.

5 Summary and conclusions

We have introduced an approach to synthetic aperture imaging with intensity-only measurements that exploits the available diversity in the illuminations. The images have the same quality as when full phase information is available. There are two stages in our approach. First, we recover field cross-correlations over pairs of frequencies at each measurement location \mbox{\boldmath{x}}_{j} using intensity-only data. Thus, the phases are recovered up to a location dependent factor $e^{\imath\theta_{j}}$, which is independent of the frequency. In the second stage, which is the main contribution of this paper, we introduce an algorithm that recovers the missing phases up to a single global factor $e^{\imath\theta_{1}}$ for all locations and all frequencies. We can then image with any method that uses full (phase and amplitude) data. We explore this approach with broadband SAR in the $50$GHz regime in an imaging setup that is used in security scanning equipments.

Appendix A Illumination strategies

We discuss here an illumination strategy for recovering the phase cross-correlations from time-resolved intensity measurements at different frequencies using a single transmit/receiver element. In this protocol at each transmit/receiver location we need to record measurements corresponding to $3S-2$ illuminations. We explain next the proposed protocol.

We want to recover the cross-correlation matrix

[TABLE]

using diversity of illuminations and the polarization identity

[TABLE]

[TABLE]

When intensities are recorded, all the quantities $|\cdot|^{2}$ in the right hand side of (21) and (22) are known.

Remarking that

[TABLE]

we deduce that we only need to compute the phase cross-correlations $m^{j}_{l1}m^{j}_{1l^{\prime}}$ which can be obtained from the polarization identity (21)-(22) provided $3S-2$ measurements. Indeed we can determine $m^{j}_{l1}$, for $l=2,\ldots,S$ using illuminations \mbox{\boldmath{e}}_{l}, \mbox{\boldmath{e}}_{l+1} and \mbox{\boldmath{e}}_{l-{\imath}1} and we also need to measure $m^{j}_{11}$.

Appendix B Measurement process

We describe here the measurement process that allows us to compute the elements $m^{j}_{ll^{\prime}}$ of the cross-correlation matrix \mbox{\boldmath{M}}^{j} corresponding to frequencies $\omega_{l}$ and $\omega_{l^{\prime}}$. When two signals of distinct frequencies $\omega_{l}$ and $\omega_{l^{\prime}}$ are sent simultaneously to probe the medium, the intensity at the receivers oscillates at the difference frequency $\omega_{l}-\omega_{l^{\prime}}$. Therefore, our imaging system should be able to resolve intensities that oscillate at frequencies of the order of the total bandwidth $B$.

Indeed, the intensity of the slowly varying terms is of the form

[TABLE]

We remark that both $I_{0}$ and $I_{\pi/2}$ oscillate with a slow frequency $\omega_{l}-\omega_{l^{\prime}}$. We assume that we can measure the integral of these terms over some time interval $\Delta t\ll\frac{1}{B}$. This leads to a linear system of the form

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

By solving system (23) subject to the constraint

[TABLE]

we determine the phases $\Delta\phi^{j}_{ll^{\prime}}$ of the elements $m^{j}_{ll^{\prime}}=|m^{j}_{ll^{\prime}}|e^{\Delta\phi^{j}_{ll^{\prime}}}$, $l\neq l^{\prime}$. Then, it is straightforward to obtain the phase differences

Acknowledgments

The work of M. Moscoso was partially supported by Spanish grant MICINN FIS2016-77892-R. The work of A.Novikov was partially supported by NSF grant NSF DMS-1813943. The work of G. Papanicolaou was partially supported by AFOSR FA9550-18-1-0519. The work of C. Tsogka was partially supported by AFOSR FA9550-17-1-0238 and AFOSR FA9550-18-1-0519.