Exact Imaging of Extended Targets Using Multistatic Interferometric Measurements

TL;DR

This paper introduces a new wave-based multistatic interferometric imaging method using Generalized Wirtinger Flow (GWF) that guarantees exact recovery of extended targets, enabling super-resolution imaging.

Contribution

It develops a deterministic lifted forward model satisfying GWF conditions for exact multistatic radar imaging, with bounds on pixel spacing and sample complexity.

Findings

GWF guarantees exact recovery under certain conditions.

The method achieves super-resolution in multistatic radar imaging.

Applicable to acoustics and geophysics imaging problems.

Abstract

In this paper, we present a novel approach that can exactly recover extended targets in wave-based multistatic interferometric imaging, based on Generalized Wirtinger Flow (GWF) theory [1]. Interferometric imaging is a generalization of phase retrieval, which arises from cross-correlation of measurements from pairs of receivers in multistatic configuration. Unlike standard Wirtinger Flow, GWF theory guarantees exact recovery for arbitrary lifted forward models that satisfy the restricted isometry property over rank-1, positive semi-definite (PSD) matrices with a sufficiently small restricted isometry constant (RIC). To this end, we design a deterministic, lifted forward model for interferometric multistatic radar satisfying the exact recovery conditions of the GWF theory. Our results quantify a lower limit on the pixel spacing and the minimal sample complexity for exact multistatic radar imaging via GWF. We provide a numerical study of our RIC and pixel spacing bounds, which shows that GWF can achieve exact recovery with super-resolution. While our primary interest lies in radar imaging, our method is also applicable to other multistatic wave-based imaging problems such as those arising in acoustics and geophysics.

Equations

$$f_i(\omega )=$$

$$\omega \in [{\omega }_c-B/2,{\omega }_c+B/2]\subset \mathbb{R}$$

$${\varphi }_i(\mathbf{x})=\left|\mathbf{x}-{\mathbf{a}}_i^r\right|+\left|\mathbf{x}-{\mathbf{a}}^t\right|$$

$$A_i(\mathbf{x},\omega )=\frac{J_i(\mathbf{x},\omega )J_t(\mathbf{x},\omega )}{\left|\mathbf{x}-{\mathbf{a}}_i^r\right|\left|\mathbf{x}-{\mathbf{a}}^t\right|}$$

$$d_{i,j}(\omega )={\int }_{D\times D}{\mathrm{e}}^{\mathrm{i}\omega /c_0{\phi }_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime })}{A}_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )\tilde{\rho }(\mathbf{x},{\mathbf{x}}^{\prime })d\mathbf{x}d{\mathbf{x}}^{\prime }$$

$${\phi }_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime })=\left|\mathbf{x}-{\mathbf{a}}_i^r\right|+\left|\mathbf{x}-{\mathbf{a}}^t\right|-\left|{\mathbf{x}}^{\prime }-{\mathbf{a}}_j^r\right|-\left|{\mathbf{x}}^{\prime }-{\mathbf{a}}^t\right|,$$

$$A_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )=A_i(\mathbf{x},\omega )A_j^{\ast }({\mathbf{x}}^{\prime },\omega )$$

$$\tilde{\rho }(\mathbf{x},{\mathbf{x}}^{\prime })=\rho (\mathbf{x}){\rho }^{\ast }({\mathbf{x}}^{\prime })$$

$${\phi }_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime })\approx \left|{\mathbf{a}}_i^r\right|-\left|{\mathbf{a}}_j^r\right|-⟨{\hat{\mathbf{a}}}_i^r,\mathbf{x}⟩+⟨{\hat{\mathbf{a}}}_j^r,{\mathbf{x}}^{\prime }⟩-⟨{\hat{\mathbf{a}}}^t,\mathbf{x}-{\mathbf{x}}^{\prime }⟩$$

$$A_{i,j}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )\approx {\alpha }_{i,j}:=\frac{C_t^2{J}_i({\mathbf{x}}_o,\omega )J_j({\mathbf{x}}_o,\omega {)}^{\ast }}{\left|{\mathbf{a}}_i^r\right|\left|{\mathbf{a}}_j^r\right|{\left|{\mathbf{a}}^t\right|}^2}$$

$$d_{i,j}({\omega }_m)=⟨{\mathbf{L}}_i^m,\mathbf{\rho }⟩⟨{\mathbf{L}}_j^m,\mathbf{\rho }{⟩}^{\ast }=\mathrm{tr}\left({\mathbf{L}}_j^m({\mathbf{L}}_i^m{)}^H\tilde{\mathbf{\rho }}\right)$$

$${\mathbf{L}}_i^m=[{\mathrm{e}}^{-\mathrm{i}{\omega }_m/c_0{\varphi }_i({\mathbf{x}}_k)}{A}_i{]}_{k=1}^K,i=1,...N.$$

$$\mathbf{d}=\frac{1}{\sqrt{M\left(\genfrac{}{}{0ex}{}{N}{2}\right)}}[{\mathbf{d}}_{1,2}^T,\dots ,{\mathbf{d}}_{N-1,N}^T{]}^T$$

$$\mathbf{d}=\mathcal{F}(\tilde{\mathbf{\rho }})$$

$$\mathbf{d}=\mathbf{F}\overline{\mathbf{\rho }}$$

$${\Phi }_{i,j}^{k,k^{\prime },l,l^{\prime }}=$$

$$+⟨{\hat{\mathbf{a}}}^t,{\mathbf{x}}_k-{\mathbf{x}}_{k^{\prime }}-{\mathbf{x}}_l+{\mathbf{x}}_{l^{\prime }}⟩.$$

$${\Vert \mathbf{d}\Vert }_2^2={\Vert \mathcal{F}\tilde{\mathbf{\rho }}\Vert }_2^2=\frac{{\sum }_{i<j}|{\alpha }_{i,j}{|}^2}{\left(\genfrac{}{}{0ex}{}{N}{2}\right)}$$

$$\times \tilde{\rho }({\mathbf{x}}_k,{\mathbf{x}}_{k^{\prime }})\tilde{\rho }({\mathbf{x}}_{l^{\prime }},{\mathbf{x}}_l))$$

$$\mathcal{K}(\Phi )=\frac{\sin \left[\left({\omega }_c^{\prime }+\frac{B}{2}\right)\frac{\Phi }{c_0}\right]-\sin \left[\left({\omega }_c^{\prime }-\frac{B}{2}\right)\frac{\Phi }{c_0}\right]}{B\frac{\Phi }{c_0}},$$

$$(1-\delta )\Vert \tilde{\mathbf{\rho }}{\Vert }_F^2\le \frac{1}{M}\Vert \mathcal{F}(\tilde{\mathbf{\rho }}){\Vert }^2\le (1+\delta )\Vert \tilde{\mathbf{\rho }}{\Vert }_F^2,$$

$$\hat{\mathbf{X}}=\frac{1}{M}{\mathcal{P}}_S\left(\mathrm{Re}\{{\mathcal{F}}^H(\mathbf{d})\right)\},$$

$$\mathcal{J}(\mathbf{\rho })=\frac{1}{2}\Vert \mathcal{F}(\mathbf{\rho }{\mathbf{\rho }}^H)-\mathbf{d}{\Vert }_2^2,$$

$$\nabla \mathcal{J}({\mathbf{\rho }}_k)=\frac{1}{M}{\mathcal{P}}_S\left(\mathrm{Re}\{{\mathcal{F}}^H({\mathbf{e}}_k)\}\right){\mathbf{\rho }}_k,$$

$$\frac{2\mu }{\alpha }\le \frac{(1-{\delta }_1)h^2({\delta }_1)}{(1+{\delta }_1)c^2({\delta }_1)},$$

$$\underset{{\omega }_c^{\prime }\to \infty }{\lim }\mathcal{K}({\Phi }_{i,j}^{k,k^{\prime },l,l^{\prime }})=\{\begin{array}{ll}0& {\Phi }_{i,j}^{k,k^{\prime },l,l^{\prime }}\ne 0\\ 1& {\Phi }_{i,j}^{k,k^{\prime },l,l^{\prime }}=0.\end{array}$$

$$\underset{{\omega }_c^{\prime }\to \infty ,S\to \infty }{\lim }$$

$$=$$

$$\times \tilde{\rho }({\mathbf{x}}_k,{\mathbf{x}}_{k^{\prime }})\tilde{\rho }({\mathbf{x}}_{l^{\prime }},{\mathbf{x}}_l)=0$$

$$M\ge \mathcal{O}\left(\frac{L}{{\Delta }_{res}}\right),$$

$${\lambda }_c=\frac{2\pi c_0}{{\omega }_c^{\prime }}$$

$$\delta \le \frac{2\pi }{A}\frac{2{\lambda }_c\sqrt{L{\Delta }_{Res}}}{{\Delta }^2\cos \varphi \sqrt{\cos \varphi }}+\mathcal{O}\left(\frac{K}{(N/A{)}^2}{\lambda }_c^{-3/2}\right)$$

$${\Delta }_{res}=2\pi \frac{c_0}{2B},\text{and}\Delta =\frac{L}{\sqrt{K}}.$$

Full text

I Introduction

I-A Motivation and Objective

In this paper, we study the exact reconstruction of complex scenes in the context of multistatic interferometric imaging. Interferometric imaging is a close relative of phaseless imaging where, in lieu of self-correlated, intensity only data, we have pairwise cross-correlated data that introduces a phase component. This work establishes Generalized Wirtinger Flow (GWF), a computationally efficient interferometric imaging method developed in [1], as a theoretical framework for exact multistatic imaging of complex scenes, while relating its recovery guarantees to the imaging system parameters. To this end, we design a deterministic and underdetermined measurement model satisfying the GWF’s sufficient condition for exact recovery. In addition, we show that it is possible to obtain exact reconstruction at resolutions smaller than Fourier-based methods and provide minimum order of measurements sufficient to guarantee such reconstruction.

The recently developed GWF algorithm is inspired by standard Wirtinger Flow (WF) [2] developed for the generalized phase retrieval problem. WF is a computationally efficient alternative to lifting based methods [3, 4]. The GWF algorithm extends the standard WF to interferometric inversion problems, and identifies a sufficient condition for exact recovery for arbitrary measurement models, characterized over the lifted domain. Hence, unlike standard WF theory which guarantees exact recovery for specific random measurement models, GWF theory guarantees exact recovery for a general class of inverse problems including random and deterministic models that abide by a single condition. In particular, the sufficient condition requires the lifted forward map to satisfy the restricted isometry property (RIP) for rank-$1$, positive semi-definite (PSD) matrices with a sufficiently small restricted isometry constant (RIC). To the best of our knowledge, our work is the first in which a deterministic and underdetermined forward model satisfying RIP for rank-$1$ PSD matrices in the lifted domain has been designed.

We provide two outcomes that unify the imaging problem with the mathematical theory of GWF. First, we determine the minimum pixel spacing to satisfy the sufficient condition for exact recovery via the GWF algorithm. Our lower bound depends on the imaging system parameters, thereby, quantifies the range of values and imaging scenarios for exact recovery guarantees to hold. For common radar imaging parameters spanning passive and active imaging modalities, this fundamental lower bound outperforms the range resolution limit of Fourier-based imaging methods for sufficiently small scenes.

Secondly, we determine the sample complexity in the order of the number of unknowns to be reconstructed. Unlike the classical results from electromagnetic scattering theory which study the degrees of freedom of scattered fields via their spatially band-limited nature and Nyquist sampling [5, 6], our analysis is based on the discrete problem, with a geometry consisting of sparsely distributed, static, terrestrial receivers. Furthermore, our sampling complexity result directly relates to the exact recovery guarantees of our algorithm through the properties of the lifted forward map, rather than the accuracy of interpolating the sampled scattered field. Notably, our result establishes that the exact recovery guarantees of GWF hold for a lifted forward model that is underdetermined. Hence, we specify a multistatic measurement model with an optimal order of measurements for interferometric wave-based imaging.

I-B Related Work and Advantages of GWF

Interferometric techniques have a rich history in acoustic, geophysical, and electromagnetic wave-based imaging. Cross correlations are frequently deployed as a fundamental formulation for passive modalities [7, 8, 9, 10, 11], in which the received ambient signal originates from a source of opportunity, such as a wireless communication signal, digital TV signal or FM radio for radar imaging. On the other hand, cross-correlations were proposed for active imaging problems to mitigate the effects of statistical fluctuations in scattering media [12, 13], clutter [14, 15, 16], and phase errors in the correlated linear transformations [17, 18, 19, 20]. More recently, inversion from phaseless scattered fields were proposed [21, 22, 23], to fully eliminate the need of coherent data acquisition in various modalities to cut implementation costs [24], or evade fundamental issues in maintaining phase coherence over long synthetic apertures [25]. Here, we motivate our approach for interferometric multi-static radar imaging via its advantages over several conventional and modern methods in the literature.

I-B1 Passive Radar

A popular method for passive imaging is the time difference of arrival (TDOA) backprojection [7, 9, 26, 27, 28, 29]. Although they are computationally efficient, TDOA backprojection is based on certain assumptions on the scatterers [7] that are not applicable for realistic scenes and can produce undesirable background artifacts [30].

To mitigate this problem, methods based on lifting have been adapted to the interferometric measurement model for passive imaging [11, 31]. These methods are inspired by the convex semi-definite programming approaches to phase retrieval [3, 4, 32, 33], which reformulate inversion as a low rank matrix recovery (LRMR) problem. Convexification has the added advantage that LRMR is known to have theoretical exact recovery guarantees under certain conditions on the lifted forward map [34]. However, these advantages come at the cost of increasing the dimension of the inverse problem, and hence introduce several limitations. Specifically, as a result of lifting, LRMR suffers from limitations on spatial sampling of the imaging grid due to high computational complexity and demanding memory requirements [11, 1].

GWF is a non-convex optimization approach that operates fully on the original signal domain, thus avoids lifting the problem at implementation and provides computational and memory efficiency over LRMR methods. Unlike TDOA/FDOA backprojection, GWF guarantees exact recovery without additional prior knowledge or limiting assumptionson the scene. Furthermore, the exact recovery guarantees afforded by LRMR [34] require more stringent conditions on the lifted forward map than that of GWF [1].

I-B2 Active Radar

For active imaging, there exists a rich literature of methods on general multistatic geometries involving distributed antennas [35, 36] or arrays [37]. These include time reversal and beamforming [38, 39], subspace methods (MUSIC [37, 40], linear sampling [41, 42, 43]), and iterative optimization schemes [44, 45, 46].

Time-reversal, beamforming and MUSIC have wide use in array imaging problems. These methods assume that the scatterers in the scene of interest are point-like and the number of measurements are greater than number of scatterers in the scene [40, 47]. This is in stark contrast to our GWF framework, in which no such assumptions are needed.

Linear sampling methods were devised to extend the applicability of subspace methods to the reconstruction of extended targets in the far field and can recover the boundaries of extended objects [41, 42]. Similar to our imaging system geometry, linear sampling methods consider a scenario that the receivers and transmitters fully encircle the scene of interest in the far field. However, the method degrades considerably when the aperture angle is less than $2\pi$ radians [48]. Our GWF result quantifies the impact of the aperture angle directly in the sufficient condition for exact recovery. Hence, GWF has applicability when aperture angle to the scene is limited.

Regularized iterative reconstruction approaches, such as total variation (TV) [44] and $\ell$-1 regularization [45, 46], have shown to achieve edge preservation. However, regularized iterative reconstruction approaches, in general, do not offer a theoretical exact recovery guarantee. Notably, TV regularization, while convex, is known to have multiple non-trivial minimizers. In addition, the TV regularizer does not have a closed form proximity operator, hence iterative reconstruction requires an inner optimization problem at each iteration. Similar problems also exist with $\ell_{1}$ regularization due to existence of a tuning parameter, which is heuristically determined. More importantly, the $\ell_{1}$ regularizer is based on strong sparsity assumptions on the scene, which is not applicable to realistic scenes. GWF, on the other hand, offers exact recovery guarantees for complex, realistic scenes with low computational complexity per iteration.

I-C Organization of the Paper

In Section II, we describe the signal model for interferometric multistatic radar. Section III presents our main results which establish how imaging system parameters of multistatic radar controls the RIC for rank-$1$ real-valued PSD matrices. Section IV describes the simulated experiments performed to verify our results in Section III. Section V concludes the paper.

II Signal Model

II-A Received Data Model

Let $N$ be the number of receivers each deployed at different spatial locations $\mathbf{a}_{i}^{r}$, $i=1,\dots,N$, where subscript $i$ denotes the $i$-th receiver and superscript $r$ denotes receiver. Assume a single transmitter located at $\mathbf{a}^{t}$. Furthermore, without loss of generality, we assume that the ground topography is flat. Thus, each spatial location in three-dimensional space is represented as $\mathbf{x}=[\bm{x},0]$ where $\bm{x}\in\mathbb{R}^{2}$. Under these assumptions and that only the scattered field is being measured under the Born approximation, the received signal at each receiver for multistatic radar can be modeled as  [49]

[TABLE]

where

[TABLE]

is the bistatic phase function, $D\subset\mathbb{R}^{2}$ is the support of the scene, $\omega$ is the fast-time frequency variable, $\omega_{c}$ is the center frequency, $B$ is the bandwidth, $c_{0}$ is the speed of light, $\rho$ is the target/scene reflectivity function; and $A_{i}$ is the amplitude function given by

[TABLE]

with $J_{i}$, and $J_{t}$ being the receiver and transmitter antenna beampatterns.

II-B Correlated Measurements

Given the data model (1), we consider the interferometric data, i.e. fast-time cross-correlation of the measurements at pairs of different receivers. Furthermore, we make the assumption that $\left|J_{t}(\mathbf{x},\omega_{m})\right|\approx C_{t}\in\mathbb{R}^{+}$. In other words, we assume that the transmitted waveform has a flat spectrum, and that the scene remains in the $-3$ dB beam-width of the illumination pattern. This is typical of radar waveforms and waveforms of opportunity such as phase shift keying (PSK) modulation found in orthogonal frequency-division multiplexing (OFDM) common among digital communications. Using (1)–(3), the correlated measurements can be modeled as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

with $(\cdot)^{*}$ denoting complex conjugation. We call $\tilde{\rho}$ the lifted version of $\rho$ or the Kronecker scene.

We next make the small-scene and far-field approximation and approximate the phase term in (5) as

[TABLE]

where $\hat{\mathbf{a}}$ is the unit vector in the direction of $\mathbf{a}$, and  (6) as

[TABLE]

where $\bm{x}_{o}$ denotes the center of the scene, with $\left|J_{i}(\bm{x}_{o},\omega)\right|\approx C_{i}$ over the observed frequency band. We assume that the support of the scene is discretized into $K$ discrete spatial points, $\{\bm{x}_{k}|\ k=1,\dots,K\}$ and define $\bm{\rho}=[\rho(\bm{x}_{1}),\dots,\rho(\bm{x}_{K})]^{T}$. We further assume that the support of $\omega$ is discretized into $M$ samples, $\Omega=\{\omega_{m}|m=1,\dots,M\}$ so that $\bm{d}_{i,j}=[d_{i,j}(\omega_{1}),\dots,d_{i,j}(\omega_{M})]^{T}$, $\omega_{m}=\omega_{c}-B/2+\frac{m-1}{M}B$.

We write (4) as

[TABLE]

where

[TABLE]

Let

[TABLE]

be the full vectorized data scaled by the number of correlated measurements. (10) shows that the data vector $\mathbf{d}$ is linear in $\tilde{\bm{\rho}}$, the Kronecker scene, while it is non-linear in $\bm{\rho}$. Thus, the data vector can be written as

[TABLE]

where $\mathcal{F}$ is a linear mapping from $\mathbb{R}^{K\times K}$ to $\mathbb{C}^{M\binom{N}{2}}$. Alternatively, if $\overline{\bm{\rho}}$ is the column-wise vectorization of $\tilde{\bm{\rho}}$,

[TABLE]

where $\mathbf{F}$ is a complex-valued matrix of size $M\binom{N}{2}\times K^{2}$, whose rows are formed by row-wise vectorization of the matrix $\mathbf{L}_{i}^{m}(\mathbf{L}_{j}^{m})^{H}$.

III Exact Multistatic Wave-Based Imaging

In this section, we are concerned with identifying the imaging system parameters, i.e., design of the measurement vectors $\mathbf{L}_{i}^{m}$, $i=1,\dots,N$, $m=1,\dots,M$, so that the lifted forward map $\mathcal{F}$ satisfies the sufficient condition proved in [1] for exact recovery via the GWF algorithm. We establish all the results presented in this section under the following assumption.

Assumption 1**.**

Let

[TABLE]

Then, we assume that $\frac{B}{2Mc_{0}}\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}\ll 2\pi$ for all $(i,j,k,k^{\prime},l,l^{\prime})$ where $B$ is the bandwidth of the received signal, $M$ is the number of frequency samples, and $c_{0}$ is the speed of light in a vacuum.

Assumption 1 is used to make small angle approximation in the proof of Lemma 1 below. This assumption implies that the number of frequency samples needed depends on the bandwidth of the transmitted waveform and the maximum value of $\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}$, which depends on the size of the scene and the placement of the receivers. As later seen in (25), this assumption is easily satisfied if the scene is sufficiently small.

We next introduce the following lemma that expresses the kernel of the operator $\mathcal{F}$ in terms of sinc functions. This lemma is used in proving Propositions 1 and 2, and for the main result in Theorem 1.

Lemma 1**.**

Suppose Assumption 1 holds. Then, the $2$-norm of the data, $\mathbf{d}$ can be written as

[TABLE]

where the phase term $\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}$ is as in (15) and

[TABLE]

with $\omega_{c}^{\prime}=\omega_{c}-\frac{B}{2M}$.

Proof.

See Appendix A-A. ∎

III-A GWF Framework for Interferometric Imaging

For establishing the GWF as an exact interferometric wave-based imaging framework, we study the restricted isometry property (RIP) of the lifted forward map $\mathcal{F}$ over the set of rank-$1$, PSD matrices.

Definition 1**.**

Let $\mathcal{F}:\mathbb{C}^{K\times K}\rightarrow\mathbb{C}^{M\binom{N}{2}}$ denote the lifted forward model provided in (13). Then $\mathcal{F}$ satisfies the restricted isometry property over rank-1, positive semi-definite matrices of the form $\tilde{\bm{\rho}}=\bm{\rho}\bm{\rho}^{H}$, with a restricted isometry constant (RIC)-$\delta$, if

[TABLE]

for any $\bm{\rho}\in\mathbb{C}^{K}$, where $\|\cdot\|_{F}$ denotes the Frobenius norm.

Notably, we consider a domain of $\bm{\rho}\in\mathbb{R}^{K}$ in optimization via GWF, for which the exact recovery guarantees of [1] directly apply. Thereby, the GWF algorithm for interferometric imaging is summarized as follows:

• •

Inputs: Interferometric measurements $d_{ij}(\omega_{m})$ and measurement vectors $\mathbf{L}_{i}^{m}$, $i=1,\cdots N$.

• •

Initialization: Backproject the interferometric measurements to the lifted domain, i.e., form an estimator for the Kronecker scene as:

[TABLE]

and keep its rank-1 approximation $\lambda_{0}\bm{\rho}_{0}\bm{\rho}_{0}^{H}$, where $\mathcal{P}_{S}$ is the projection operator onto the set of symmetric matrices. The initialization step consists of the outer product of the two measurement vectors for each of the $M$ samples, resulting in $\mathcal{O}(MK^{2})$ multiplications, followed by an eigenvalue decomposition with $\mathcal{O}(K^{3})$ complexity.

• •

Iterations: Initializing with $\bm{\rho}_{0}$, perform gradient descent updates as $\bm{\rho}_{k+1}=\bm{\rho}_{k}-\frac{\mu_{k}}{\|\bm{\rho}_{0}\|^{2}}\nabla\mathcal{J}(\bm{\rho}_{k})$, having

[TABLE]

which yields

[TABLE]

where $(\mathbf{e}_{k})_{ij}^{m}=(\left(\mathbf{L}_{i}^{m}\right)^{H}\bm{\rho}_{k}\bm{\rho}_{k}^{H}\mathbf{L}_{j}^{m}-d_{ij}(\omega_{m})),\ \forall i\neq j$, $m=1,\cdots,M$.

Each iteration requires the following operations:

Computing and storing the linear terms $(\mathbf{L}_{i,j}^{m})^{H}\bm{\rho}_{k}$, requiring $M$ number of $K$ multiplications for each, resulting in $\mathcal{O}(MK)$ multiplications. 2. 2.

Computing the error by cross correlating linear terms, requiring $\mathcal{O}(M)$ multiplications. 3. 3.

Multiplication of the linear terms $(\mathbf{L}_{i,j}^{m})^{H}\bm{\rho}_{k}$ and the error $e_{ij}^{m}$ for each $m=1,\cdots M$, requiring $\mathcal{O}(M)$ multiplications. 4. 4.

Multiplication of the result in $3$ with vectors $\{\mathbf{L}_{i}^{m}\}_{m=1}^{M}$ and $\{\mathbf{L}_{j}^{m}\}_{m=1}^{M}$, requiring $\mathcal{O}(MK)$ multiplications.

These operations result in $\mathcal{O}(MK)$ multiplications for each iteration.

In particular, [1] establishes that exact recovery of a ground truth $\bm{\rho}_{t}\in\mathbb{C}^{K}$ is guaranteed upto a global phase factor, if the forward operator for the lifted Kronecker scene, satisfies the RIP over the set of rank-$1$, PSD matrices (i.e., $\tilde{\bm{\rho}}$) with RIC of less than $0.214$. Furthermore, starting from the initial estimate computed from  (19), gradient descent iterations minimizing  (20) using  (21) converges geometrically to the true solution at a rate $1-\kappa$, where $\kappa$ is upper bounded by

[TABLE]

where $h\leq c$ are positive constants solely depending on $\delta$ [1].

III-B Asymptotic Result

As a stepping stone for our main result, we begin by showing the asymptotic isometry of $\mathcal{F}$ defined in (13), as $\omega_{c}\rightarrow\infty$ and $N\rightarrow\infty$. Following our asymptotic analysis of the kernel of $\mathcal{F}$, we characterize its RIP over rank-1, PSD matrices in the non-asymptotic regime. As a result of our non-asymptotic analysis, we derive an upper bound on the restricted isometry constant that is controlled by the imaging system parameters.

Despite its limited use in practice, our initial asymptotic result offers a valuable benchmark for the non-asymptotic case. Notably, it justifies assessing how the isometry of $\mathcal{F}$ is perturbed over the set of rank-1, PSD matrices when the central frequency $\omega_{c}$, and the number of receivers $N$ are finite. We specifically make use of this perspective in establishing our main result, by analytically evaluating elements $\bm{\rho}\bm{\rho}^{H}$ in the range of $\mathcal{F}^{H}\mathcal{F}$. Furthermore, it characterizes the expected limiting behavior of our upper bound estimate on the RIC-$\delta$.

The following proposition shows that in the asymptotic regime, i.e., as $\omega_{c}$ gets large, $\mathcal{F}$ becomes a delta function with respect to the phase term $\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}$.

Proposition 1**.**

Under Assumption 1, we have

[TABLE]

Proof.

See Appendix A-B. ∎

Given Proposition 1, the next proposition shows that in the limit as $\omega_{c}\rightarrow\infty$ and $N\rightarrow\infty$, $\mathcal{F}$ is an isometry.

Proposition 2** (Asymptotic Isometry of $\mathcal{F}$ for large $\omega_{c}$ and $N$).**

Under Assumption 1, we have

[TABLE]

Proof.

See Appendix A-C. ∎

Since in the asymptotic regime $\mathcal{F}$ is an isometry, we can deduce that the RIC over rank-$1$, PSD should become small as $\omega_{c}$ and $N$ get large. This motivates us to find an upper bound on the rank-$1$, PSD RIC constant in the non-asymptotic regime in terms of the imaging parameters. In the next subsection, we establish this upper bound.

III-C Non-asymptotic Result

Before we introduce our main theorem, we introduce two further assumptions.

Assumption 2**.**

The scene is enclosed by a square with side $L$ and sampled regularly on a square grid. The coordinate system is centered at the middle of the square. Hence, $\bm{x}=[x_{1},x_{2}]^{T}\in[-L/2,L/2]\times[-L/2,L/2]$ with $\sqrt{K}$ samples in both $x_{1}$- and $x_{2}$-axis and $L=\sqrt{K}\Delta$ where $\Delta$ is the pixel spacing.

Under Assumption 2, it is easy to see that the phase term $|\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}|$ is upper bounded by $4L\sqrt{2}$ for any selection of $i,j,k,k^{\prime},l,l^{\prime}$.

Then, for Assumption 1, letting $\Delta_{res}=2\pi\frac{c_{0}}{2B}$ be the range resolution given by the Fourier-based methods the small angle approximation holds to high accuracy if

[TABLE]

since $\max_{i,j,k,k^{\prime},l,l^{\prime}}|\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}|=\mathcal{O}\left(L\right)$. For instance, $M\geq 5.8\frac{L}{\Delta_{res}}$ corresponds to a $<1\%$ error for the $\mathrm{sinc}$ approximations in Lemma 1.

Assumption 3**.**

The receivers are isotropic and lie on a circular arc equidistant from each other and to the center of the coordinate system. Let $A\in(0,2\pi]$ be the aperture of the multistatic system. Then, the azimuth angles of the look-directions are multiples of $A/N$. 2. 2.

*All receivers and the transmitter are located at the same height. Let $\phi$ be the elevation angle in radians. Then, $\hat{\mathbf{a}}^{r}_{i}=[\cos\phi\cos\theta_{i},\cos\phi\sin\theta_{i},\sin\phi]$ where $\theta_{i}=\frac{Ai}{N}$, $i=0,\dots,N-1$ are the azimuth angles of the receivers’ look-directions. * 3. 3.

The transmitter is located on the $x_{1}$-axis. Hence, $\hat{\mathbf{a}}_{t}=[\cos\phi,0,\sin\phi]^{T}$.

Assumption 3 allows us to make integral approximation to a Riemann sum in the proof of Theorem 1 (see Appendix A-D). The approximation error is then incorporated into the result of Theorem 1. Note that the assumption on the location of the transmitter is not essential, but is there for convenience.

We now state our non-asymptotic result in the following theorem, which establishes an upper bound on the rank-$1$, PSD RIC for the data model presented in (13), in terms of the underlying imaging parameters.

Theorem 1** (RIC of the Lifted Forward Mapping of Multistatic Imaging).**

Let

[TABLE]

be the wavelength corresponding to the center frequency. Then, under Assumptions 2, 3, and Lemma 1, we have the following upper bound on the restricted isometry constant $\delta$ of $\mathcal{F}$ over rank-1, PSD matrices:

[TABLE]

where the order is a small constant and

[TABLE]

Proof.

See Appendix A-D. ∎

As provided in  (22) and explained in [1], $\delta$ directly controls the convergence rate of GWF iterates. As such, bound in (27) establishes that the convergence behavior of GWF for multistatic imaging depends on system parameters such as the center frequency $\omega_{c}$, the bandwidth $B$, the number of receivers $N$, the number of unknowns $K$, as well as the side length $L$ of the scene.

Remark 1*.*

Observe that $N$ has a higher order than $\lambda_{c}^{-1}$ in the second term in (27). Hence, our RIC upper bound estimate tends to [math] as $\omega_{c}\rightarrow\infty$, $N\rightarrow\infty$, consistent with our asymptotic isometry result for $\mathcal{F}$. Specifically, the first term in (27) captures the perturbation from the limit when the central frequency is finite, whereas the second term characterizes the perturbation due to having finite number of receivers. In fact, the second term directly arises from the closed form error of a Riemann sum approximation to an integration over look directions of the receivers.

Remark 2*.*

The Riemann sum error behaves in an inverted manner to the first term with respect to the central frequency of the transmitted signal is increased, given a fixed imaging aperture and number of look directions. This is indeed an expected outcome, as the data collection manifold corresponds to a larger area of the 2D Fourier spectrum of the scene as the central frequency is increased while the imaging aperture is fixed. As a result, $N$ number of look directions corresponds to a poorer discretization of the data collection manifold, and the factor of $\lambda_{c}^{-3/2}$ in the second term in (27) relates directly to this phenomenon.

Using the decoupled nature of our upper bound estimate on the RIC, we quantify the minimal pixel spacing at which the exact recovery guarantees of GWF can hold.

Corollary 1** (Resolution Bound).**

Suppose we have sufficiently many receivers, i.e., $N^{2}\gg K$, such that the second term in (27) is negligible. Then GWF guarantees exact recovery if

[TABLE]

Proof.

Assuming $N^{2}\gg K$, the second term in (27) in the upper bound of $\delta$ vanishes. Recall that exact recovery is guaranteed via GWF if $\delta$ is less than or equal to $0.214$. Upper bounding the RIC bound in (27), we have

[TABLE]

The rest follows by rearranging (30). ∎

Notably, even with $N\rightarrow\infty$, (29) is the absolute best resolution at which exact multi-static imaging is possible by GWF. Hence, Corollary 1 yields a fundamental bound for the pixel spacing in designing realizable imaging systems with finite number of receivers.

The resolution bound of Corollary 1 corresponds to the super-resolution regime when reconstructing small scenes in both active, and passive scenarios, as depicted in Figure 1(a) and Figure 1(b), respectively. Note that as $L$ gets large, the lower bound eventually becomes greater than the range resolution limit of the Fourier-based methods. This is in agreement with our theoretical arguments, which are established under a small scene approximation. It should also be stressed that our lower bound abides by the sufficient condition for exact recovery, but it is not a necessary one. Therefore, while recovery of scenes at a higher resolution than $\Delta_{res}$ may still be possible via GWF, it is not covered by the theory in [1].

Additionally, the sufficient number of receivers for (29) to hold is $N^{2}\geq\mathcal{O}(K)$. Since $M=\mathcal{O}(L)$ by (25), this implies that super-resolution imaging via GWF requires a sample complexity of at least $MN^{2}=\mathcal{O}(K^{3/2})$. We reduce this complexity result by the following corollary, which quantifies the minimal sample requirement for exact multi-static imaging via GWF at a fixed pixel spacing that abides the lower bound of Corollary 1 .

Corollary 2** (Sample Complexity).**

Given the final result of Theorem 1, exact multistatic imaging condition for GWF is satisfied at the following sample complexity:

[TABLE]

Proof.

Reorganizing the upper bound on $\delta$ in Theorem 1, we have

[TABLE]

where $c_{1}$, $c_{2}$ are $\mathcal{O}(1)$ as functions of $K$ and $N$. Now, for any fixed pixel spacing $\Delta$, we have $L=\mathcal{O}(\sqrt{K})$. Thus,

[TABLE]

for some $\tilde{c}_{1},\tilde{c}_{2}=\mathcal{O}(1)$. Observe that $K^{1/4}$ factor in the first term of the left-hand side of (33) is non-vanishing and hence at best yields the RIC upper bound of $\tilde{\delta}=\mathcal{O}(K^{1/4})$. Now from Assumption 1, we have $M=\mathcal{O}(L)$. Thus, the minimal sample complexity in which the RIC upper bound is in the order of $K^{1/4}$ is achieved when $N^{2}=\mathcal{O}({K}^{3/4})$. Therefore,

[TABLE]

when $MN^{2}=\mathcal{O}(K^{5/4})$. ∎

In addition to the minimal sample complexity,  Corollary 2 yields a rate at which the algorithm performance deteriorates. Clearly, from (34), our ability to fine sample the scene while attaining the exact recovery guarantees of GWF for multi-static imaging depends on the dimension of the problem, at a rate $K^{1/4}$, or equivalently, $\sqrt{L}$. This, again, is consistent with our theoretical arguments as we derive our results through a small scene approximation.

The fact that the upper bound of $\delta$ has a non-vanishing $K^{1/4}$ factor reveals an interesting phenomenon that is also observed in the performance of spectral initialization in phase retrieval literature, even when the measurement vectors are random. This degradation with the increasing dimension of the unknown is not captured in the probabilistic analysis with random measurement vectors, yet is indeed a significant issue which forms the basis for sample truncation in computing the initialization and gradient estimates [50].

Specifically for deterministic, wave-based multistatic imaging problems,  Corollary 2 necessitates a system design such that the controllable constants in (34) sufficiently suppress the $K^{1/4}$ factor. This promotes GWF as a highly applicable method in passive imaging scenarios where the range resolution is limited, or in active imaging scenarios where small, isolated extended targets are being imaged, with possible extensions and applications in spot-light mode synthetic aperture radar [51].

IV Numerical Simulations

In this section, we provide several numerical simulations demonstrating veracity of the theory presented in Section III. The following multistatic set-up is common to all simulations presented in this section, and conforms to the assumptions laid out in Section III.

There is a single transmitter located at $[15.8,0,0.25]$ km. 2. 2.

The transmitted waveform has unit amplitude frequency spectrum. 3. 3.

Varying number of receivers are distributed equidistant on an arc of a circle of radius $10$ km from the scene center at a height of $0.25$ km. 4. 4.

The scene of interest is square with flat topography.

Figure 2 illustrates the multistatic set-up used in this section. Note that the illustration is not to scale.

The figure-of-merit we use throughout is the mean square error (MSE) of the reconstructed scene. This is computed by taking the per pixel difference between the true scene and the reconstructed scene and averaging the squares of the differences.

In all experiments presented, our data is synthetically generated using our received signal model in  (1) under the single-scattering assumption. In Sections IV-A, IV-B, and IV-C, a single parameter is varied in each set of experiment while all other relevant parameters are fixed. The parameters are chosen in the active and passive imaging ranges. Figure 3 shows the scene used for all experiments in the subsequent sections. Finally, in Section IV-D, we provide simulations that depict the performance of GWF in non-ideal conditions, namely under additive noise at low signal-to-noise ratios (SNR), and discretization mismatch between the ground truth, and the reconstructed image.

IV-A Effect of Number of Receivers

The first series of numerical experiments are designed to verify the effect of the number of receivers on the performance of GWF reconstruction. In (27), the second term involves the square of the number of receivers, $N^{2}$, in the denominator. Thus, we expect the number of receivers to have significant effect on the quality of the reconstruction. To verify the effect of the number of receivers on the reconstruction, we ran a series of simulations with varying number of receivers while fixing all other relevant parameters in active or passive radar regimes.

Figure 4 shows the MSE of the resulting reconstruction versus the number of receivers for active and passive imaging. Blue solid line is the result for the active case while black dashed line is for the passive case. For the active case, the bandwidth was held at $B=50$ MHz with the center frequency at $\omega_{c}=10$ GHz for Fourier-based range resolution of $\Delta_{res}=3$ m. For the passive case, $B=10$ MHz and $\omega_{c}=1.9$ GHz for $\Delta_{res}=15$ m. The pixel spacing was chosen such that it was smaller than the Fourier-based range resolution for each case. Namely, $\Delta=2.4$ m and $\Delta=12$ m for the active and passive cases, respectively. The number of unknowns was held constant at $K=625$ for both cases. The GWF algorithm was performed for $4000$ iterations for comparison purposes. Since the RIC directly affects the rate of convergence of GWF, we expect to see smaller MSE as the number of receivers grows. This behavior is clearly present in both the active and passive cases as can be readily observed in Figure 4. In both cases, we observed exact convergence behavior from $10$ receivers onward. However, as expected, the convergence rate is slower with smaller number of receivers.

As a visual confirmation of the experimental verification, sample reconstructions at two different number of receivers ($12$ and $24$) is provided in Figures 5 and 6 for active and passive regimes, respectively.

IV-B Effect of Bandwidth/Range Resolution

Next we examine the effect of the bandwidth on the convergence behavior. Both terms in (27) includes square root of $\Delta_{res}$, the range resolution, in the numerator. This suggests that there is an inverse relationship between the bandwidth and RIC. Similar to above, we test the effect of bandwidth on the convergence behavior of GWF algorithm, and we ran a series of GWF reconstruction on the same scene while varying the bandwidth and holding other relevant parameters fixed. The number of receivers used for the experiments was fixed at $N=18$. All other parameters were held to the same values as in the previous subsection.

Figure 7 summarizes the result of these experiments. Figure 7(a) shows the bandwidth vs. MSE curve for active case. We varied the bandwidth in the range of $30$ MHz to $70$ MHz. Figure 7(b) shows the same curve for the passive case where the bandwidth was varied between $6$ MHz and $24$ MHz. Examining the two figures, we clearly see that higher bandwidth results in smaller MSE, and hence faster convergence to exact solution. This agrees with the theoretical bound in (27). As before, we provide visual confirmation in form of sample reconstructions in Figures 8 and 9 for active and passive regimes, respectively.

IV-C Effect of Center Frequency

The first term of (27) is inversely proportional to the center frequency of the transmitted waveform and as such we expect the center frequency to improve the convergence behavior of GWF as the center frequency gets larger. We examined numerically, the effect of center frequency on the exact reconstruction and the convergence rate by, again, running a series of numerical simulations where we varied the center frequency while keeping other relevant variables constant. To minimize the effect of the second term on the RIC and better evaluate the impact of central frequency in the super-resolution regime, we increased the number of receivers used in these experiments to $N=32$ for both cases.

Figures 10(a) and 10(b) show the results of simulated experiments for active and passive scenarios, respectively. For active case, we varied the center frequency in the range between $0.5$ GHz and $15$ GHz. For the passive case, the range was restricted to $0.1$ GHz to $3$ GHz to reflect realistic values for sources of opportunity. In both cases, we observe a behavior of downward trend in MSE as the center frequency increases, albeit, not as drastic as in other parameters. This is due to the fact that the central frequency appears inverted in the two terms of the RIC upper bound. The decaying trend of our experiments agrees with the notion that the order constant of the second term in the RIC upper bound adequately suppresses $\lambda_{c}^{-3/2}$, as $N^{2}=\mathcal{O}(K)$ proves to be sufficient for super-resolution.

Notice, however, that in the active case, larger center frequency value is needed to achieve similar performance as in the passive case. This is attributable to the fact that the first term is proportional to $\sqrt{L\Delta_{res}}/\Delta^{2}$. With the active parameters, this term is approximately $8$ times that of the passive case. Thus, the center frequency needs to be higher to compensate for the difference. Figures 11 and 12 show sample reconstructions at two different center frequencies for active and passive regimes, respectively.

IV-D Effect of Additive Noise

Next, we evaluate the robustness of the proposed method for interferometric multi-static radar imaging, in passive, and active cases. As before, we use a flat spectrum signal with $10$ MHz bandwidth centered around $1.9$ GHz frequency in the passive, and a $50$ MHz bandwidth centered around $10$ GHz frequency in the active experiments. We incorporate additive, zero-mean white Gaussian noise on the linear signal model at the receivers, and consider SNR levels varying from $-20$ to $30$ dB.  Figure 13 demonstrates the reconstruction MSE with respect to the received signal SNR with $5$ dB increments, with errors averaged over $10$ realizations. The MSE curves indicate that GWF is robust to additive noise at the receivers in both active and passive cases, and the performance of the algorithm degrades predictably with decreasing SNRs, with steady decay in MSE as the SNR at the receivers improve.

It should be noted that the correlation operation amplifies the noise variance, hence the interferometric measurements processed in the experiments have lower SNR than the specified levels at the receivers. Nonetheless, we observe that the reconstruction performance of GWF degrades gracefully as SNR decrases, as  Figure 14,  Figure 15 and  Figure 16 demonstrate that GWF is capable of producing highly accurate imagery with [math] dB in the received signals at the antennas. However with SNRs below $-10$ dB the algorithm performance degrades noticeably, which motivates filtering or sample truncation at the implementation of GWF in low SNR scenarios.

V Conclusion

In this paper, we utilize GWF theory developed in [1] for exact multistatic imaging of extended targets by designing the underlying imaging parameters such that the sufficient condition for exact recovery is satisfied. Our work has two significant contributions.

1. Unlike the state-of-the-art interferometric inversion methods based on LRMR, GWF avoids lifting the problem. As a result, it is computationally efficient and does not require large memory allocations, making it suitable for practical applications.
2. We demonstrate that the underlying imaging parameters can be designed so that the RIP over rank-1, PSD matrices is satisfied by a deterministic lifted forward model.

We first show the asymptotic isometry of the lifted forward model, $\mathcal{F}$, of interferometric multistatic radar, as the center frequency and the number of receivers go to infinity. We then proceed with estimating the deviation from the asymptotic behavior when imaging parameters are finite, and derive an upper bound for the RIC of $\mathcal{F}$ over the set of rank-1, PSD matrices. Hence, we identify the relation of imaging parameters to the sufficient condition of exact recovery. Using the RIC upper bound, we determine a lower limit for pixel spacing to achieve exact recovery. This limit is superior to the Fourier-based range resolution for sufficiently small scenes. Furthermore, we determine the minimal sample complexity needed for RIC upper bound to be sufficiently small, hence identify the practical requirements for reconstruction when designing a multistatic imaging system. In our numerical simulations, we evaluate the impact of the imaging parameters in our upper bound estimate of RIC in reconstruction performance, verify our theoretical results, and finally assess the performance of GWF to additive noise at low SNRs.

For future work, we will study the robustness of our method with respect to deviations from our imaging setup, such as non-equi-distant locations, or non-circular configurations of receivers. In addition, we will investigate extensions of our theory to the case involving additive noise and outliers, moving target imaging, as well as implementation of our method using real scattering data.

Appendix A Proofs

A-A Proof of Lemma 1

We first examine the $2$-norm of the data. For a rank-$1$ $\tilde{\bm{\rho}}$, we have that $\left\|\tilde{\bm{\rho}}\right\|_{F}^{2}=\left\|\bm{\rho}\right\|_{2}^{4}$. We can also rewrite

[TABLE]

Thus, from (11), (8), and (9) we have

[TABLE]

Similarly, we have that

[TABLE]

where $\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}$ is as in (15).

Then, under Assumption 1, we have that

[TABLE]

where $\omega_{c}^{\prime}=\omega_{c}-\frac{B}{2M}$. The second line is from geometric sum and the last line is from small angle approximation.

Using (38) and changing the order of sum, and denoting $\left|\alpha_{i,j}\right|^{2}=\frac{\left|C_{i}\right|^{2}\left|C_{j}\right|^{2}}{\left|\mathbf{a}_{i}^{r}\right|^{2}\left|\mathbf{a}_{j}^{r}\right|^{2}\left|\mathbf{a}^{t}\right|^{4}}$, we have

[TABLE]

We can split (39) into two parts as

[TABLE]

Having a real-valued $\tilde{\bm{\rho}}$, we rewrite the latter term in (40) as

[TABLE]

We can further rewrite $\mathcal{K}(\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}})$ (41) using trigonometric identity as

[TABLE]

which proves the claim.

A-B Proof of Proposition 1

First we express $\mathcal{K}$ as

[TABLE]

where

[TABLE]

Given (43), it suffices to prove that

[TABLE]

This can be proved using similar machinery to proving the delta function limit for sequence of scaled sinc functions. We prove the first equation (45). The second equation follows similarly.

Let $f\in S(\mathbb{R})$ be a smooth test function, where $S(\mathbb{R})$ is the Schwartz space. Then we need to prove that

[TABLE]

where

[TABLE]

Let $\epsilon>0$ and break-up the integral into two parts.

[TABLE]

The first part of the integral is

[TABLE]

We compute the first part of (50). Integrating by parts,

[TABLE]

Since $f\in S(\mathbb{R})$, taking the limit as $\omega_{c}^{\prime}\rightarrow\infty$, (51) goes to [math]. Similarly, we can see that the second integral in (50) also goes to zero. Thus,

[TABLE]

Now, the second integral in (49) can be rewritten as

[TABLE]

The first integral can be rewritten as

[TABLE]

Since $f\in S(\mathbb{R})$, we can use Riemann-Lebesgue lemma to conclude that (54) goes to zero as $\omega_{c}^{\prime}\rightarrow\infty$ [52].

For second integral, we have

[TABLE]

After change of variables $u=\left(\omega_{c}^{\prime}+\frac{B}{2}\right)\frac{x}{c_{0}}$, we have

[TABLE]

as $\omega_{c}^{\prime}\rightarrow\infty$. We can use similar argument for where

[TABLE]

A-C Proof of Proposition 2

Without loss of generality, let the receivers and transmitters have common elevation angle $\phi$ such that

[TABLE]

where $\theta_{i}$ is the azimuth angle of the $i$-th receivers look-direction, $\theta_{t}$ is the azimuth angle of the transmitter look-direction and $\phi$ is the elevation angle. Furthermore, we have that for any $k$ and $k^{\prime}$

[TABLE]

where $\theta_{k,k^{\prime}}$ is the angle of the vector $\bm{x}_{k}-\bm{x}_{k^{\prime}}$. Then we have that

[TABLE]

Thus, for the non-diagonal terms where $k\neq k^{\prime}$, $l\neq l^{\prime}$ we have that $\Phi_{i,j}^{k,k^{\prime},l,l^{\prime}}=0$ if

[TABLE]

For fixed $k,k^{\prime},l,l^{\prime}$ and $i$, there are at most $2$ values of $\theta_{j}$’s for which (62) is satisfied. Furthermore, we know that $\left|\alpha_{i,j}\right|^{2}$’s must be bounded. Thus, by Proposition 1, for each fixed $k,k^{\prime},l,l^{\prime}$ where $k\neq k^{\prime}$ and $l\neq l^{\prime}$ we have that

[TABLE]

Now taking the limit as $N\rightarrow\infty$, we have the desired result.

A-D Proof of Theorem 1

We want to upper bound the following

[TABLE]

where, $\left|\alpha_{i,j}\right|^{2}\approx\alpha=\frac{\left|C\right|^{4}}{(\bm{a}_{1}^{r})^{4}(\bm{a}^{t})^{2}}$ by Assumption 3. Without loss of generality we set $\alpha=1$. We begin by noting that

[TABLE]

where

[TABLE]

Thus, fixing $l$, $l^{\prime}$, and $k$, we have convolution between $\mathcal{G}$ and $\rho$. Let

[TABLE]

We take the Fourier Transform of $\mathcal{G}_{i,j}$ and $\rho$ to represent the convolution. Denoting, $\hat{\mathcal{G}}_{i,j}$ as the Fourier Transform of $\mathcal{G}_{i,j}$, we have

[TABLE]

To compute $\hat{\mathcal{G}}_{i,j}$, we first rewrite $\mathcal{G}_{i,j}$ as

[TABLE]

Let $\bm{x}_{k}=[x_{1}^{k},x_{2}^{k}]^{T}$, $\bm{\omega}=[\omega_{1},\omega_{2}]^{T}$ and $\theta_{i}$ be the azimuth angle of the $i$-th receiver’s look-direction. Then, given (69), the Fourier Transform of $\mathcal{G}$ and using the assumption that $\hat{\bm{a}}_{t}=[1,0]^{T}$,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Noting that $R$ is only non-zero where $\frac{(\omega_{c}^{\prime}-B/2)}{c_{0}}\cos\phi(\cos\theta_{i}+1)\leq\omega_{1}\leq\frac{(\omega_{c}^{\prime}+B/2)}{c_{0}}\cos\phi(\cos\theta_{i}+1)$, and $\omega_{c}^{\prime}\gg B/2$, we approximate (71) as

[TABLE]

Next, we note that

[TABLE]

Thus, interchanging the sum and the integral in (68), and plugging in (70) we have

[TABLE]

where

[TABLE]

Now, by employing Cauchy-Schwartz, we have

[TABLE]

where

[TABLE]

By Jensen’s inequality, (80) becomes

[TABLE]

Noting that for any fixed $\omega_{1}$,

[TABLE]

we have

[TABLE]

We use Jensen’s inequality once more to get

[TABLE]

Next, approximating the sum over $\theta_{j}$ as an integral, we have

[TABLE]

where $E_{R}$ denotes the Riemann sum error, $A$ is the aperture of look directions in the imaging setup.

We first consider the inner integration over $\omega_{1}^{\prime}$, where $\bm{\omega}^{\prime}=[\omega_{1}^{\prime},\theta_{j}]$. Using Cauchy-Schwartz and Jensen’s inequalities,

[TABLE]

Computing the first integral in (88), we get

[TABLE]

Consider the second integral in (88). Since the integrand is strictly positive, from the $\theta_{j}$ integration we have

[TABLE]

We now make the following change of variables

[TABLE]

Computing the Jacobian, we get

[TABLE]

Thus, setting $\bm{\omega}^{\prime\prime}=[\omega_{1}^{\prime\prime},\omega_{2}^{\prime\prime}]$, the upper bound in (90) becomes

[TABLE]

where the last identity follows from Parseval’s theorem. Hence, we obtain the upper bound on our integral approximation, i.e., the first term in (87) as

[TABLE]

We next evaluate the error term of the integral approximation in (93). Namely, the midpoint Riemann-sum approximation error is upper bounded as

[TABLE]

where $Q$ is defined as

[TABLE]

using the definition of $\tilde{R}(\omega_{1}^{\prime})$.

Since the integrand is the squared absolute value of the Fourier transform of the reflectivity function evaluated at frequencies $\bm{\omega}^{\prime}$, (96) can equivalently be written as

[TABLE]

which is the form of a Frobenius inner product, where

[TABLE]

Observe that $Q$ has an upper bound that only depends on the $\ell_{2}$ norm of the underlying scene as

[TABLE]

and that $\underset{\theta}{\mathrm{max}}\|\frac{\partial^{2}}{\partial\theta^{2}}\tilde{\mathbf{f}}\|_{F}$ has only dependence on the imaging system, hence yields a universal upper bound for any scene reflectivity function in the specific imaging geometry. Evaluating the integral in (98) and having $\omega_{c}\gg B$, $\tilde{f}$ can be approximated under the narrow-band assumption as

[TABLE]

Denoting $z=\frac{\omega_{c}}{c_{0}}\gamma(\theta)\cdot(\bm{x}_{l}-\bm{x}_{l^{\prime}})$, the second order derivative of each component in $\tilde{f}$ with respect to $\theta$ is obtained as

[TABLE]

where

[TABLE]

[TABLE]

By definition, the amplitude of sinc-function derivatives in (101) have a decay rate of $1/z$, whereas $z$ and its derivatives in (101) grow with an order of the scene dimension $L$, as $\mathcal{O}(\frac{\omega_{c}}{c_{0}}L)$. Hence, evaluating the squared integral for the Frobenius norm, the growth is at a maximal order of $\mathcal{O}\left((\frac{\omega_{c}}{c_{0}})^{3}L^{3}\right)$. We thereby obtain a final approximate upper bound on the universal constant $Q$ as

[TABLE]

Putting the terms derived in (103) and  (94) into (86) and (77), we have

[TABLE]

for the definition of $\delta$ in Definition 1, using the fact that $\|\bm{\rho}\|^{4}=\|\tilde{\bm{\rho}}\|_{F}^{2}$ for the rank-1, PSD element $\tilde{\bm{\rho}}=\bm{\rho}\bm{\rho}^{H}$.

Re-organizing the terms using (26) and (28) defined in the statement of Theorem 1, and the fact that

[TABLE]

we finally obtain

[TABLE]

which completes the proof.