Performance Analysis of Convex LRMR based Passive SAR Imaging

# Performance Analysis of Convex LRMR based Passive SAR Imaging

Eric Mason and Birsen Yazıcı, Yazıcı and Mason are with the Department of Electrical, Computer and Systems Engineering, Rensselaer Polytechnic Institute, 110 8th Street, Troy, NY 12180 USA E-mail: B. Y: yazici@ecse.rpi.edu, Phone: (518)-276 2905, Fax: (518)-276 6261. Corresponding author.This work was supported by the Air Force Office of Scientific Research (AFOSR) under the agreement FA9550-16-1-0234, and by industrial and government membership fees to the Center for Surveillance Research, a National Science Foundation Industry and University Cooperative Research Program (I/UCRC) under agreement IIP-1539961.
###### Abstract

Passive synthetic aperture radar (SAR) uses existing signals of opportunity such as communication and broadcasting signals. In our prior work, we have developed a low-rank matrix recovery (LRMR) method that can reconstruct scenes with extended and densely distributed point targets, overcoming shortcomings of conventional methods [1]. The approach is based on correlating two sets of bistatic measurements, which results in a linear mapping of the tensor product of the scene reflectivity with itself. Recognizing this tensor product as a rank-one positive semi-definite (PSD) operator, we pose passive SAR image reconstruction as a LRMR problem with convex relaxation.

In this paper, we present a performance analysis of the convex LRMR-based passive SAR image reconstruction method. We use the restricted isometry property (RIP) and show that exact reconstruction is guaranteed under the condition that the pixel spacing or resolution satisfies a certain lower bound. We show that for sufficiently large center frequencies, our method provides superior resolution than that of Fourier based methods, making it a super-resolution technique. Additionally, we show that phaseless imaging is a special case of our passive SAR imaging method. We present extensive numerical simulation to validate our analysis.

## I Introduction

A passive synthetic aperture radar (SAR) uses signals of opportunity such as communication or broadcasting signals to form an image. Passive radars offer several advantages including reduced cost, simplicity of hardware, increased stealth and efficient use of electromagnetic spectrum. With the proliferation of transmitters of opportunity, there is a growing interest in passive radar applications [2, 3, 4, 5, 6, 7, 8, 9, 10].

Widely used existing methods for the passive SAR problem can be categorized into two classes: passive coherent localization (PCL) [11, 12, 13, 7, 14, 15] and time-difference-of-arrival (TDOA) based backprojection [2, 16, 3, 5, 4, 17, 18]. These approaches rely on certain assumptions and imaging configurations that limit their applicability. Specifically, PCL requires a direct line-of-sight to a transmitter of opportunity and knowledge of transmitter location, and TDOA-based backprojection relies on a point target assumption on the scene. In [1, 19], we presented a low-rank matrix recovery (LRMR)-based alternative method that overcomes these shortcomings. This method can image scenes with extended or densely distributed point targets, and does not require direct-line-of-sight to the transmitter or knowledge of the transmitted waveform [1].

We assume that two or more spatially separated moving receivers acquire scattered field data from a scene of interest illuminated by transmitters of opportunity. We next cross-correlate received signals from different receivers. The correlation process induces a linear mapping from the tensor product of the scene reflectivity with itself to cross-correlated measurements. We linearize the problem by using the well-known \qqlifting approach, by noting that the tensor product of the scene reflectivity with itself is the kernel of a rank-one positive semi-definite (PSD) operator [20, 21, 22, 23]. Thus, we approach image reconstruction as a rank minimization problem and pose it using the well-known nuclear-norm as a convex relaxation [24, 25].

In this paper, we present performance analysis of convex LRMR-based passive SAR imaging method. Specifically, we show that under the condition that the distance between scatterers satisfies a certain lower bound, the recovery is exact, i.e., nuclear-norm minimization algorithm converges and recovers the location and reflectivities of scatterers exactly. The condition on the distance between the scatterers determines the highest achievable resolution for exact recovery. Given the bandwidth and operating frequencies of typical illuminators of opportunity, the achievable resolution is an order of magnitude higher than that of Fourier based techniques, making the convex LRMR-based passive SAR imaging a super-resolution technique.

Our analysis is based on the use of the restricted isometry property (RIP), one of several tools in LRMR theory used to establish exact reconstruction guarantees when solving nuclear-norm minimization problems [24, 26]. This property quantifies how close the forward model is to an isometry in terms of the restricted isometry constant (RIC), when the forward map is restricted to a set of rank-r matrices. We show that under certain conditions the RIC associated with the passive SAR forward model is small enough to guarantee that the solution to the nuclear-norm minimization problem is exact. We use the method of stationary phase in determining the RIC.

Nuclear-norm minimization problem can be solved in a variety of ways [24, 27, 28, 29]. Here, we pose it as a trace-minimization, PSD constrained problem with an affine data-fidelity constraint, and use Uzawa’s method to solve the resulting saddle point problem111 This approach is different from the trace regularized least-squares formulation we used in [1]. The new formulation allows us to establish its equivalence to the nuclear norm-minimization problem that the RIP based recovery results apply.. We then prove convergence to a unique minimizer, which recovers the exact unknown scene reflectivity when the RIP holds.

Additionally, we show that phaseless imaging is a special case of our LRMR-based passive imaging method. Assuming the two receivers are co-located our passive SAR forward model reduces to an auto-correlation. This type of configuration is known as phaseless imaging [20, 21, 30, 31]. Recovery results have been obtained for active phaseless imaging using a rectangular array in [23]. Our analysis covers the phaseless imaging case and generalizes the result in [23] to arbitrary imaging geometries and passive imaging. It also provides a sufficient condition that proves it to be a super-resolution technique. Additionally, we show that cross-correlation results in a smaller RIC than that of auto-correlated measurements, thus requiring lower center frequencies to obtain exact reconstruction.

The rest of this paper is organized as follows: in Section II we present the forward model and discretization of the problem. In Section III, we formulate a set of optimization problems for passive SAR imaging , present an optimization algorithm and use the relevant LRMR theory used in our analysis. Section IV presents our main result and a sufficient condition for exact recovery. Section V provides an extension of our analysis to phaseless imaging. In Section VI we discuss the implications of our analysis and the validity of our assumptions. Section VII presents numerical simulations. Section VIII concludes the paper.

## Ii Passive Imaging Forward Model

We assume there are two airborne receivers flying over a scene of interest, illuminated by a transmitter of opportunity located at . We denote the receiver trajectories by where represents the slow-time variable. An example of this configuration is displayed in Figure 1.

be the fast-time temporal frequency variable, where is the bandwidth and is the center frequency. is the baseband frequency set. The location on the surface of the earth is given by , where and is a known ground topography. denotes the ground reflectivity. Table I lists all major notation used throughout the paper.

Under the start-stop and Born approximations, assuming waves propagate in free space, the fast-time temporal Fourier transform of the received signal at the -th receiver due to a stationary transmitter located at can be modeled as [32]

 fi(ω,s)≈Li[~ρ](ω,s):=∫e−iωϕi(s,\bix,y)/c0×Ai(ω,s,\bix)~ρ(\bix)d\bix,  i=1,2, (1)

where

 ϕi(s,\bix,y)=|y−x|+|x−γi(s)|, (2)

and

 Ai(ω,s,\bix)=Ji(ω,\bix)(4π)2|x−γi(s)||y−x|. (3)

In (3), models the antenna beam patterns and transmitted waveform.

### Ii-B Model for Correlated Measurements

We correlate the bistatic measurements acquired by two spatially separated receivers, given in (1), as follows:

 d12=f1(ω,s)f∗2(ω,s) (4)

where denotes complex conjugation. Substituting (1) into (4), assuming the range of transmitter is much larger than that of the scene (i.e. ), we model the correlated measurements as follows [1]:

 d12(ω,s)≈F12[ρ](ω,s):=∫e−iωϕ12(s,\bix,\bix′)/c0×A12(ω,s,\bix,\bix′)ρ(\bix,\bix′)d\bixd\bix′ (5)

where

 ρ(\bix,\bix′)=~ρ(\bix)~ρ∗(\bix′) (6)

and

 ϕ12(s,\bix,\bix′)=|x−γ1(s)|−|x′−γ2(s)|+^y⋅(x−x′), (7)
 A12(ω,s,\bix,x′)=J1(ω,\bix)J∗2(ω,\bix′)(4π)2|x−γ1(s)||x′−γ2(s)||y|2. (8)

In (7), denotes the unit vector in the direction .

We define where denotes the tensor product. Clearly, is a rank-one positive semi-definite (PSD) operator with kernel . We refer to the kernel of as the Kronecker scene. Since is rank-one, the scene reflectivity is the only eigenfunction of , and can be recovered exactly from . While is non-linear in scene reflectivity , it is linear in . Therefore, our goal is to recover the Kronecker scene using .

In the rest of this paper, for notational brevity, we drop the subscripts from and . Furthermore, we represent the kernel of as

 F(ω,s,\bix,\bix′)=e−iωϕ(s,\bix,\bix′)/c0A(ω,s,\bix,\bix′). (9)

### Ii-C Discretized Model

We discretize the scene reflectivity using the standard pixel basis and write

 d(ω,s)≈N∑k=1N∑k′=1F(ω,s,\bixk,\bixk′)ρ(\bixk,\bixk′) (10)

where are the entries of an rank-one positive semi-definite matrix, defined as . We vectorize this matrix by stacking the columns on top of each other, forming the vector . The data is discretized by sampling uniformly in slow-time and fast-time to form a matrix with entries . We then vectorize the data matrix by stacking the columns (each slow-time sample) in a vector, defined as . The forward model is discretized using the same sampling schemes for the scene and data, this results in a matrix . The columns and rows are ordered. In particular, for the -th row, the -th entry is given as

 [Fmp]kk′=e−iωmϕ12(sp,\bixk,\bixk′,^y)/c0A12(ωm,sp,\bixk,\bixk′). (11)

Thus, we write

 ¯d=F¯ρ. (12)

We use this discrete data model in applying the LRMR theory, and for numerical implementation.

## Iii LRMR based Image Reconstruction

### Iii-a Optimization based Imaging

Since the unknown discretized Kronecker scene is a rank-one matrix we can pose image formation as the following rank minimization problem:

 P1: minimizeρ J1(ρ)=%rank(ρ)s.t.F¯ρ=¯d. (13)

This problem is non-convex and NP-hard. The difficulty of this problem has motivated development of alternative heuristics, two of which apply to PSD matrices and are referred to as the log-determinant and trace heuristics [33, 25, 34]. More recently, within the context of phaseless imaging, a penalty with a higher degree of regularization that promotes a rank-one solution has been developed [22].

While some of these alternative penalties may perform better as a result of the structure of our problem, we use the nuclear norm heuristic due to its theoretical connection to Problem [24]. This penalty is applicable to general rectangular matrices and leads to the following convex problem [24, 34]:

 P2: minimizeρ J2(ρ)=∥ρ∥∗s.t.F¯ρ=¯d, (14)

where is the nuclear-norm222Also known as the Schatten-1 norm or trace norm., defined as

 ∥ρ∥∗=∑iσi(ρ), (15)

with being the singular values of . From (14), it is clear that the convex relaxation of the rank minimization problem is the -norm of the singular values, which draws a direct analogy between cardinality minimization and rank minimization [26].

For relatively small scale problems (matrices of size less than ), the nuclear-norm minimization problem in (14) can be reformulated as a semi-definite program and solved efficiently using interior point methods [24]. However, for large scale problems, semi-definite programming fails and alternative algorithms must be used. Since the nuclear norm is non-smooth, a popular and efficient algorithmic approach is to use forward-backward splitting methods [28].

In our problem, since the unknown is positive semi-definite, an equivalent formulation of is

 P3: minimizeρ J3(ρ)=%Tr[ρ]s.t.F¯ρ=¯d,  ρ⪰0 (16)

where denotes that is positive semi-definite, and denotes the trace operator.

To solve we introduce Lagrange multipliers and use a saddle point method. However, first we incorporate the PSD constraint to the objective function using the following indicator function

 iC+(ρ)={0 ρ∈C+∞ρ∉C+ (17)

where is the set of PSD matrices. Thus, becomes

 P4: minimizeρ J4(ρ)=%Tr[ρ]+iC+(ρ)s.t.F¯ρ=¯d. (18)

To obtain a simple iterative scheme, it is useful to enforce that the solution have minimum norm. Thus, we modify as

 P5: minimizeρ J5(ρ)=%Tr[ρ]+iC+(ρ)+1λ∥ρ∥2Fs.t.F¯ρ=¯d. (19)

Next, we introduce a Lagrange multiplier to incorporate the equality constraint into the objective function and obtain

 L(ρ,ξ)=Tr[ρ]+iC+(ρ)+1λ∥ρ∥2F+⟨ξ,F¯ρ−¯d⟩. (20)

While is no longer the same as , its solution is equal to that of as . This follows from Theorem 3.1 in [27]. In our numerical simulations, we show that the solution to (20) converges to the true solution for a sufficiently large .

Under strong duality assumption, the solution to (19) can be obtained by solving the dual problem, using the Uzawa iterative procedure [27]: first update as the minimizer of (20) for fixed , then update the Lagrange multiplier with a gradient ascent update. The final iterative scheme then becomes:

 {ρk  =PC+(Rs(FHξk)−λI)ξk+1=ξk+βk(F¯ρk−¯d) (21)

where is a projection onto the PSD cone, is the step size, is an identity matrix, , , and .

The convergence of (21) to the solution of (18) is stated below in Theorem 1.

###### Theorem 1.

Suppose that the sequence of step sizes obeys , where is the Lipshitz constant333Note that the can be calculated as the largest singular value of . of . Then the sequence generated by (21) converges to the unique solution of , defined in (18).

###### Proof.

See Appendix B. ∎

### Iii-B Recovery Guarantees

We now provide the LRMR theory that states the conditions under which the rank-one solution to is unique and the same as the solution to [24]. These recovery guarantees are developed under the assumption that the measurement matrix or forward model satisfies the RIP generalized to low-rank matrices given below [24].

###### Definition 1.

Restricted Isometry Property Let be a linear operator. Assume, without loss of generality, . For every , define the -restricted isometry constant to be the smallest depending on , such that

 (1−δr)∥X∥F≤∥A[X]∥2≤(1+δr)∥X∥F (22)

holds for all matrices , of rank of at most , where is the Frobenius norm.

The restricted isometry constant (RIC), , can be viewed as a measure of how close a linear operator is to an isometry when restricted to a domain of matrices with rank at most . The relevant theorems generally state that when the RIC is less than unity there is a unique solution to the rank minimization problem, , and for an appropriately smaller this solution is the same as that of .

We begin our review of LRMR theory by considering the data model where . Two recovery results, based on the assumption that the underlying forward model satisfyies the RIP condition with small enough RIC, are given in Theorems 2 and 3 [24].

###### Theorem 2.

Let and suppose for some integer . Then, is the only matrix of rank at most satisfying the matrix equation , and hence the only solution of .

###### Proof.

See proof of Theorem 3.2 in [24]. ∎

###### Theorem 3.

Suppose and . Then, is the unique solution to .

###### Proof.

See proof of Theorem 3.3 in [24]. ∎

Theorem 2 states that for , the solution to in (13) is unique and equal to the exact solution . Theorem 3 states that for , is the unique and exact solution to in (14). Since when , Theorem 3 implies Theorem 2. From this, we conclude that when satisfies RIP with , the solutions of and are the same and equal to the true Kronecker scene .

In [23] the RIC condition in Theorem 3 is improved and exact recovery is shown for a rank-one matrix when . We summarize this result in the following Theorem.

###### Theorem 4.

Let with . Assume that . Then is the unique solution to .

###### Proof.

See proof of Theorem 3.4 in [23]. ∎

## Iv Performance Guarantees for LRMR based Passive SAR Imaging

In this section, we prove an asymptotic recovery result for LRMR-based passive SAR. In summary, our result shows that exact reconstruction is possible if the target spacing is larger than a fraction of the wavelength corresponding to center frequency. Lemma 6 states the precise conditions required for exact recovery. Theorem 7 then shows that the RIC is small enough to satisfy the conditions stated in Theorems 2, 3 and 4, proving that the solution to is unique and the same as the solution to .

 ∥Fρ∥22= M∑m=1P∑p=1∣∣ ∣∣N∑k=1N∑k′=1F(ωm,sp,\bixk,\bixk′)ρ(\bixk,\bixk′)∣∣ ∣∣2 (23) = N2∑k,k′=1N2∑l,l′=1(M∑m=1P∑p=1F(ωm,sp,\bixk,\bixk′) ×F∗(ωm,sp,\bixl,\bixl′))ρ(\bixk,\bixk′)ρ∗(\bixl,\bixl′).

In order to determine RIC, we seek an estimate of the inner summation,

 K(¯\bixk,¯\bixl):=M∑m=1P∑p=1F12(ω,s,¯\bixk)F∗12(ω,s,¯\bixl) (24)

where we define and .

Assuming sufficiently many samples in fast-time frequency and slow-time, we use (5) to approximate (24), as follows444In (25), we apply a change of variables and evaluate the integral in baseband .:

 K(¯\bixk,¯\bixl)≈∫S×Ωe−i(ωc+ω)θ(s,¯\bixk,¯\bixl)/c0B(ω,s,¯\bixk,¯\bixl)dsdω (25)

where

 θ(s,¯\bixk,¯\bixl)= |xk−γ1(s)|−|xk′−γ2(s)|−|xl−γ1(s)| (26) +|xl′−γ2(s)|+^y⋅(xk′−xk−xl′+xl)

and

 B(ω,s,¯\bixk,¯\bixl)=1(4π)4|y|4B1(s,¯\bixk,¯\bixl)B2(ω,¯\bixk,¯\bixl) (27)

with

 B1(s,¯\bixk,¯\bixl)= 1|xk−γ1(s)||xk′−γ2(s)| (28) ×1|xl−γ1(s)||xl′−γ2(s)|,

and

 B2(ω,¯\bixk,¯\bixl)=J1(ω,\bixk)J∗2(ω,\bixk′)J∗1(ω,\bixl)J2(ω,\bixl′). (29)

In order to obtain RIC, we need to estimate the integral in (25). To do so, we evaluate (25) for those where is zero for all and for those where is non-zero for almost all using the method of stationary phase. However, before we can estimate the integral, it is necessary to obtain a lower bound on . This lower bound requires a minimum achievable resolution stated in Lemma 5.

We define the set of all possible pixel locations as . We then divide this set into two disjoint sets and . Specifically, these two sets satisfy and . Additionally, is zero for all indices in and for all . For each index in , is zero for at most two values of 555This can be seen by setting over and solving for the roots of the resulting second degree polynomial.. Hence, is non-zero over for almost all .

By studying the behaviour of (25) restricted to each of these sets, we obtain the results in Lemmas 5 and 6.

###### Lemma 5.

Assume the ground topography is flat666The result can be extended to non-flat topography. However, for simplicity of argument only flat topography case is considered. and . Let denote the minimum spacing between two scatterers, and be the angle between and its orthogonal projection onto the ground plane. Assume

 Δx≫c0ωc|cosαy−1|. (30)

Then, in (26) satisfies

 min(k,k′,l,l′)∈I2{ωcc0∣∣θ(s,¯\bixk,¯\bixl)∣∣}≫1. (31)
###### Proof.

See Appendix C-A. ∎

###### Lemma 6.

Suppose

 CJ=∫ΩB2(ω,¯\bixk,¯\bixk)dω<+∞ (32)

and . Let

 W(s,¯\bixk,¯\bixl)=∫Ωe−iωc0θ(s,¯\bixk,¯\bixl)B2(ω,¯\bixk,¯\bixl)dω. (33)

Assume satisfies

 sups,¯\bixk,¯\bixl∈R ∣∣∂α¯\bixk∂β¯\bixl∂ηsW(s,¯\bixk,¯\bixl)∣∣≤Cs(1+s2)(2−|η|)/2 (34)

where is a compact subset of and is a constant that depends on . Furthermore, assume the conditions of Lemma 5 hold and the ranges of the receivers to the scene are much larger than the size of the scene . Then, as ,

 K(¯\bixk,¯\bixl)≈⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩CJΔs(4π)4|y|4L2g1L2g2(k,k′,l,l′)∈I1,1(4π)4|y|4L2g1L2g2{√2πωc×CK+O(1ω3/2c)}(k,k′,l,l′)∈I2, (35)

where is the range of the receiver from the center of the synthetic apertures and

 CK= ∑{s0|˙θ(s0,¯\bixk,¯\bixl)=0}eiωcθ(s0,¯\bixk,¯\bixl)eiπ/4¨θ(¯\bixk,¯\bixl,s0)√|¨θ(s0,¯\bixk,¯\bixl)| (36) ×W(s0,¯\bixk,¯\bixl),

with satisfying and .

###### Proof.

See Appendix C-B. ∎

These two lemmas are necessary to obtain exact reconstruction. The result of Lemma 5 is required to ensure that we are in the asymptotic regime allowing us to use the method of stationary phase in proving Lemma 6. The result of Lemma 5, in return, requires a minimum target spacing, yielding a lower bound on the achievable resolution for exact reconstruction.

Under the assumptions of Lemmas 5 and 6, we can now obtain the following performance guarantee.

###### Theorem 7.

Suppose the conditions of Lemma 5 and 6 hold. Then as , the exact locations and amplitudes of scatterers can be recovered by solving in (14).

###### Proof.

See Appendix C-C. ∎

## V Phaseless Passive Imaging and Exact Recovery

When the two receivers are collocated, i.e., the passive SAR forward model in (5) reduces to that of phaseless imaging. In this case, the phase function in (26) becomes

 θ(s,¯\bixk,¯\bixl)= |xk−γ(s)|−|xk′−γ(s)|−|xl−γ(s)| (37) +|xl′−γ(s)|+^y⋅(xk′−xk−xl′+xl).

We show that the uniqueness results obtained in [23] also follow from our analysis, generalizing it to arbitrary imaging geometries.

In the previous section, we decompose the pixel locations into two sets and so that (26) is zero on and non-zero on . In the case of phaseless imaging, these sets are different. More specifically, there is a set for which the phase function (37) is zero, but (26) is not. We define this set as , and let , and . We carry out the same analysis in Section IV using the sets and . However, for phaseless imaging it is useful to consider and separately, instead of .

As in the previous case we can obtain recovery guarantees, however, the RIC is larger when using auto-correlated measurements than cross-correlated measurements. This is a result of the fact that in the additional set in phaseless imaging

Before we state the recovery results for phaseless imaging, we briefly comment on the proofs of Lemmas 5 and 6 with and replaced with and , respectively. First, the proof of Lemma 5 remains the same, since in this case we only need to consider the difference between the distance of two different scatterer locations to the receivers. This is the case, because the only way three of the four magnitude terms in the phase cannot be equal is for one them to be equal to the fourth, and this would make the phase function zero. Thus, the approach to the proof and the resulting final condition follow from the same procedure as the proof of Lemma 5. In the case of Lemma 6, the proof remains exactly the same, despite changing the sets. Therefore, we can use these two lemmas to now state the main recovery result for phaseless imaging.

###### Theorem 8.

Suppose for all . Let and . Suppose the conditions in Lemma 5 and 6 hold with and replaced with and , respectively. Then as , the exact locations and amplitudes of the scatterers can be recovered by solving .

###### Proof.

See Appendix C-D. ∎

## Vi Discussion

We begin our discussion with the connection of the algorithmic approach in (21) to the recovery results stated in Theorems 7 and 8. The recovery results state that when the measurement map satisfies RIP with sufficiently small RIC, the solutions of and are unique and equal. This implies that the exact Kronecker scene can be recovered. While can be solved efficiently using interior point methods, this approach is not suitable for the large scale problem of SAR imaging. Therefore, in order to solve efficiently we propose the use of , which is strictly convex and has unique solution. For sufficiently large , this solution is the minimum Frobenius-norm solution of which has unique solution for sufficiently large . Therefore, when and are sufficiently large, the iteration in (21) converges to the true Kronecker scene and thus recovers the scene reflectivity exactly.

Along with provable exact reconstruction, our analysis provides a condition on the minimum resolvable target spacing. The condition is stated in Lemma 5 and suggests that convex LRMR based methods have \qqsuper-resolution capability, in the sense that it offers superior resolution than that of Fourier-based imaging. While the resolution condition depends on and , which are not user controlled parameters in passive SAR imaging, for typical terrestrial illuminators of opportunity the values of these parameters are such that super-resolution can be achieved. This is an appealing attribute of our method as resolution does not depend on bandwidth as in the case of Fourier-based methods. For example, DVB-T and WiMAX waveforms have bandwidth of - MHz and operate at around MHz and GHz center frequencies, respectively. If the transmitter elevation angle is degrees, the resolution lower bound is m for DVB-T and even smaller for WiMAX. This is far superior than the -m of resolution achievable by Fourier-based imaging.

By comparing the results of Theorems 7 and 8, it is clear that using cross-correlated measurements results in a smaller RIC than using auto-correlated measurements. This implies that required for exact reconstruction is smaller in the case of cross-correlated measurements. Smaller RIC is due to the fact that , i.e., the phase of the forward map is zero on a larger set of pixels in phaseless imaging. Hence, we conclude that a smaller RIC results from retaining phase information in the cross-correlated data. While the recovery guarantees are more easily obtained by cross-correlating measurements, auto-correlating offers incoherent imaging capability. This is useful when there is synchronization error, a common problem in radar imaging using distributed apertures [7].

We conclude this section with a discussion of the validity of the assumptions required in proving the recovery results. The main assumptions are made in evaluation of (25), stated in Lemma 6. The required assumptions apply to the imaging geometry, transmitted waveform, and antenna beam patterns. With respect to the imaging geometry, we require that the ranges of the antennas are much smaller than the size of the scene and synthetic aperture. This assumption is satisfied by all typical airborne and space-borne imaging geometries [35]. With respect to assumptions (32) and (33), they are satisfied when the antennas are sufficiently broadband and the transmitted waveform has flat spectrum or constant modulus. These are likely to hold as most antennas are designed to be broadband, and many illuminators of opportunity transmit waveforms with flat spectrum, such as DVB-T and WiMAX signals.

## Vii Numerical Simulations

In this section we validate the recovery result presented in the previous sections.

### Vii-a Experimental Set-up

We consider a scene with flat topography that is discretized into pixels. This results in a Kronecker scene of pixels, each pixel having of resolution. The scene consists of a single extended target with three different reflectivity values as shown in Figure 2. The corresponding Kronecker scene is shown in Figure 2.

We assume the scene is illuminated by a stationary transmitter located at . For flat topography, this results in a transmitter elevation angle degrees. We use system parameters corresponding to common illuminators of opportunity. Specifically, we assume a transmit pulse with a flat spectrum and MHz bandwidth, and center frequencies of MHz and GHz, corresponding to DVB-T and WiMAX signals, respectively [8, 6]. We consider a semi-circular aperture, i.e., , where the second receiver traverses the same trajectory with a constant offset of , i.e., . An illustration of this imaging configuration is displayed in Figure 3. Based on this configuration, the resolution lower bound in (30) is m and m for the MHz and GHz centrer frequencies, respectively. Thus, the m pixel spacing used in our experiments satisfies the resolution condition.

We assume isotropic transmit and receive antennas, and for simplicity assume that geometric spreading has been compensated for during data collection. This results in the amplitude function of the forward model in (8) being set to . The data at each receiver are simulated using (1), then correlated according to (4).

### Vii-B Experimental Results

#### Vii-B1 Performance Metrics and Evaluation Criteria

We denote the true and reconstructed Kronecker scene as and , and evaluate the reconstruction performance to validate exact recovery () using three error metrics:

 Ed(ρ)=∥d−F¯ρ∥2∥d∥FEρ(ρ)=∥ρ−ρ⋆∥F∥ρ∥FE~ρ(~ρ)=∥~ρ−~ρ⋆∥F∥~ρ∥F (38)

where stands for the Frobenius norm. These metrics evaluate how well the data fidelity constraint is satisfied, and how close the Kronecker scene and scene reflectivity are to the actual values. We consider success for values smaller than error with correct trace and rank. The error metrics show that we have obtained a feasible solution, and exact trace and rank confirm that we have reached a minimizer of , which is equivalent to that of . These conditions also verify that realistic center frequencies satisfy the assumptions of the recovery Theorems.

We run simulations for both cross-correlated and auto-correlated data. The auto-correlated data corresponds to passive phaseless imaging. We run the algorithm stated in (21) for iterations with to recover . While the algorithm may not have completely converged in all cases, we consider a longer running time to be impractical. We then reconstruct the scene reflectivity from by keeping the eigenvector corresponding to the largest eigenvalue.

#### Vii-B2 Reconstruction Results

In Table II we tabulate the values of the three error metrics in (38), as well as the trace and the rank of the reconstructed Kronecker scene. The true trace of the Kronecker scene displayed in 2 is . We find that using cross-correlated measurements results in exact reconstruction for both center frequencies and the reconstructed images visually match the Kronecker scene and reflectivity phantom shown in Figure 2. The rank and trace match that of the true scene, and all error metrics in (38) are below the threshold. From this we conclude that the optimal point of - has been recovered, validating our analysis.

While we obtain exact recovery for cross-correlated measurements, the algorithm fails to reconstruct the exact solution from auto-correlated or phaseless measurements. This suggests that the center frequency is not large enough for Theorem 4 to hold. Unlike our numerical results, [36, 23] demonstrate exact reconstruction. This is explained by the fact that they are using GHz or GHz center frequencies, which are not realistic operating frequencies for passive radar.

In Figure 4 we display the reconstructed images for GHz center frequency. Figures 4 and 4 display the reconstructed Kronecker scene using cross-correlated and auto-correlated measurements, respectively. The scene reflectivity corresponding to these images are shown in Figures 4 and 4, respectively. Visually, the reconstruction is exact when using cross-correlated measurements. In the case of phaseless passive imaging, it is clear that the algorithm fails to converge to the correct solution. The resulting image is significantly blurred and the true target shape and reflectivity values are not recovered.

#### Vii-B3 Algorithm Behaviour

To demonstrate the convergence behavior of our method, we plot the recovered trace trace, rank, and verses iteration number in Figure 5, for each experiment in Table II. The plot of is shown in Figure 5. We observe that in every case the error decreases monotonically, and appears to be converging to zero. When using phaseless measurements it appears that the algorithm has not fully converged, confirmed by the plots of trace in 5. Also, based on the plot of trace and rank, it appears that the algorithm is converging to an incorrect solution. This is explained by the fact that the center frequency is not large enough for the result of Theorem 8 to hold. This means the solution of is not unique-and-equal to that of . This is further verified by observing that the rank is converging to a value larger than one, which implies that the minimum norm solution of is not rank-one.

Theorem 1 states that the sequence of Kronecker scene estimates converge to the true solution, implying that and trace converge to the correct value, as observed when using cross-correlated measurements. Although, an interesting behavior is that the rank and trace increase monotonically while converging. While this may seem odd at first, it is explained by the form of the algorithm. Specifically, (21) does not explicitly solve by minimizing the trace of , but instead solves the dual problem. The solution of the dual problem lower bounds the Lagrangian at the optimal point, and since strong duality holds equality is obtained. Therefore, at each iteration the trace of the iterate increases until the lower bound converges to the optimal value of the solution.

As stated before, at the GHz center frequency Theorem 8 does not hold. Therefore, the solution of and are not the same as that of . This implies that the solution to and may not be rank-one, in fact, the solution likely has rank greater than one. Hence, the rank increases per iteration since the algorithm is converging to a solution with rank greater than one. This is further verified by the fact that when cross-correlated measurements are used the result of Theorem 7 applies and the unique solution is rank-one. Thus, at each iteration the rank of the iterate remains at one.

This illustrates the difference in the behavior of the algorithm for cross- and auto-correlated measurements. It also offers a heuristical method to determine early on if the analysis holds given the imaging parameters. If the rank increases above one, it can be assumed that the center frequency is not large enough and the algorithm can be terminated to save computational resources.

## Viii Conclusion

In this work, we presented a performance analysis of the LRMR based convex passive SAR imaging method. Using the restricted isometry property we show that exact recovery is possible for sufficiently large center frequencies. In evaluating the restricted isometry constant we determine a sufficient condition for exact recovery. This condition arises from the application of the method of stationary phase. The condition states that the distance between scatterers must be larger than a certain bound that depends on the center frequency of the transmitted waveform and the elevation angle of the transmitter. For typical transmitters of opportunity operating in the MHz to GHz range, this condition translates into an achievable resolution that is an order of magnitude higher than that of Fourier based methods used in PCL and TDOA-based backprojection making convex LRMR-based passive SAR imaging a super-resolution technique.

Furthermore, we show that our analysis applies to phaseless passive imaging, in which measurements are auto-correlated instead of cross-correlated. Our recovery results apply to phaseless imaging, however, with larger RIC and larger center frequency than those required by passive SAR with cross-correlated measurements.

To solve the imaging problem we present a projected gradient descent type algorithm and prove its convergence. We verify our analysis with extensive numerical simulations. We show that for typical operating frequencies used in passive radar, using cross-correlated measurements the scene reflectivity can be recovered exactly.

While we consider primarily delay based cross-correlated measurements, our methodology can be extended to study recovery results for Doppler-based correlations [2, 3]. Additionally, our results apply to passive synthetic aperture imaging acoustics and geophysics as well as passive imaging with sufficiently large distributed receivers [37, 38, 10].

## Appendix A Method of Stationary Phase

The stationary phase theorem states that if u is a (possibly complex-valued) smooth function of compact support in , and is a real-valued function with only non-degenerate critical points, then as

 ∫eiλφ(\bix)u(\bix)dn\bix=∑{\bix0|Dφ(\bix0)=0}(2πλ)n/2×eiλφ(\bix0)eiπ/4sgnD2φ(\bix0)√|detD2φ(\bix0)|u(\bix0)+O(λ−n/2−1). (39)

Here denotes the gradient of , and denotes the Hessian, and sgn denotes the signature of a matrix, i.e. the number of positive eigenvalues minus the negative ones. Furthermore, a point is called a non-degenerate critical point if