THE RELATIONSHIP BETWEEN BACKPROJECTION AND BEST LINEAR UNBIASED ESTIMATION IN SYNTHETIC-APERTURE RADAR IMAGING

. In this paper we investigate the relationship between two diﬀerent techniques typically used in imaging and estimation problems. We focus on synthetic-aperture radar imaging and compare the methods of backprojection (standard for imaging) and best linear unbiased estimation (BLUE). We aim to reconstruct or estimate the reﬂectivity function of an object present in a scene of interest. We ﬁnd that the estimate of the reﬂectivity (calculated using BLUE) and the reconstructed image (calculated using ﬁltered backprojection) are the same when we utilize a criterion from microlocal analysis to deﬁne the optimal backprojection ﬁlter and assume the measured data is corrupted by zero-mean independently identically distributed (white) noise. In particular we show that the microlocal criterion for the optimal backprojection ﬁlter is equivalent to the unbiased constraint present in the BLUE technique.

1. Introduction. In this paper we investigate the relationship between two ostensibly different statistical methods: backprojection in synthetic aperture imaging versus best linear unbiased estimation, a core method in classical statistical science. The best linear unbiased estimator (BLUE) coming from statistical estimation theory [7,9,8] is used to estimate parameters which summarize collected data properties. The BLUE is constrained to be an unbiased linear combination of the data. We choose the 'best' estimator using some criterion of optimality, commonly one aims to find the estimator with minimum variance or minimum mean square error.
In radar imaging we aim to estimate the reflectivity function of an object of interest which is present in a scene viewed by the radar. The radar emits electromagnetic pulses and collects measurements of the field scattered off the object and its surroundings. We use these collected scattered field measurements to estimate the unknown reflectivity function. One may think of the estimator as the reconstructed image of the scene of interest.
We relate BLUE to the method of backprojection reconstruction of the reflectivity function of an object of interest from synthetic-aperture radar (SAR) data. Backprojection is the technique of analytically inverting the collected data to find the reflectivity function. It is assumed in this case that the data can be written as an operator (namely a Fourier integral operator) applied to the reflectivity function. The inversion, or backprojection, operation aims to apply the adjoint of this 'forward' operator in order to obtain a reconstruction, or image, of the scene's reflectivity. Typically we are restricted to an approximate adjoint because we are limited by the radar's finite bandwidth and antenna aperture. We will focus on a form of filtered backprojection that utilizes microlocal analysis to define a filter that gives us an ideal image where the singularities of the object of interest appear in the correct location and with the correct orientation in the reconstructed image. We argue that in this case the image we obtain is precisely the BLUE in the case of minimizing the variance of the estimator.
Connections between the analytic technique of backprojection and techniques from estimation and detection theory has been previously investigated. The author previously looked at the comparison between the generalized likelihood ratio test (GLRT) and backprojection in [13,14]. It was found that the test statistic calculated to detect and estimate the location of a single target from data corrupted by independent identically distributed Gaussian noise using the GLRT was equivalent to the backprojection image of the target calculated using a matched filter under certain conditions. In [5] it was found that the GLRT was equivalent to taking the inner product of the Wigner distribution of the data with the Wigner distribution of the quantity one wishes to estimate. The Wigner distribution is a second order quantity which is nonlinear and clearly different from backprojection. In addition in [10] the GLRT used for estimating the location of an object from x-ray data was compared to backprojection and it was found under certain conditions that the likelihood function can be interpreted as a form of backprojection but it may not be thought of as a true image reconstruction because the kernel or filter is not chosen as it is in image reconstruction algorithms. In addition the author notes that the GLRT is a nonlinear technique which is different from backprojection based image reconstruction, a linear estimation technique. Also in [1] the general problem of detecting and localizing a point target in noise from time-harmonic scalar waves emitted from and received on a sensor array was considered. It was found that reverse-time migration (which is essentially backprojection) of the data provided the most powerful test for a given false alarm rate, using the Neyman-Pearson lemma [7], i.e. it is equivalent to the likelihood ratio test.
In this work it is found the the BLUE is more naturally related to backprojection imaging in that both are inherently linear techniques, while the GLRT is only linear under specific assumptions. In addition the filter that arises from the BLUE is equivalent to that of the backprojection filter when a criterion from microlocal analysis is used to determine the optimal filter. In [16] a similar connection was noted between the filter chosen in backprojection and the bias of the image. Also the simplifying assumptions in [13,14] are not needed to prove the relationship (for example, we do not need to assume a single target). We therefore conclude that the BLUE is the real equivalent technique to backprojection stemming from estimation and detection theory. The purpose of this work is to provide insight into the backprojection methods of SAR data analysis by drawing a relationship between this method and a central statistical method which is linear unbiased estimation.
The rest of the paper is organized as follows. In section 2 we derive the expression for SAR data from Maxwell's equations. In section 3 we discuss the method of filtered backprojection and its relationship to microlocal analysis and Fourier Integral Operator (FIO) theory. Here we also derive the optimal backprojection filter in terms of a microlocal criterion. In section 4 we derive the best linear unbiased estimator of the reflectivity function and show the equivalence to the backprojected image from section 3. We conclude with some remarks on the implications of our finding and possible future areas of study.
2. Data model. In SAR, the waves are electromagnetic and therefore their propagation is described by Maxwell's equations. We follow the data derivation present in [3,13]. In the time domain Maxwell's equations have the form where E is the electric field, B is the magnetic induction field, D is the electric displacement field, H is the magnetic intensity or magnetic field, ρ is the charge density, and J is the current density. Since most of the propagation takes place in dry air we assume that the electromagnetic properties of free space hold. Therefore we have that ρ = 0 and J = 0. We also have the free space constitutive relations which are expressed as follows: where µ 0 and 0 are the magnetic permeability and electric permittivity in free space respectively. Note we make the assumption that the scatterers are isotropic and hence 0 and µ 0 are scalars instead of tensors. This ultimately allows us to decouple the equations into three scalar equations in terms of Cartesian coordinates of R 3 .
Using these assumptions we take the curl of (1) and substitute the result into (2). Then using the triple-product identity (also known as the BAC-CAB identity) on the left-hand side we obtain Note that ∇ · E = 0 because of the free space assumption, which implies we have This is equivalent to stating that in Cartesian coordinates each component of the vector E satisfies the scalar wave equation. Using similar steps we may say that H also satisfies the wave equation in free space. In particular these expressions hold for a constant wave speed, that is c 0 = (µ 0 0 ) −1/2 . When a scatterer is present the wave speed is no longer constant and it will change when the wave interacts with the object. We may think of this as a perturbation in the wave speed, which we write as where T is known as the scalar reflectivity function. Using this variable wave speed we say that E satisfies the following wave equation: This model is often used for radar scattering despite the fact that it is not entirely accurate. The reflectivity function T really represents a measure of the reflectivity for the polarization measured by the antenna.
In order to incorporate the antennas we consider the total electric field E tot = E in + E sc which is composed of the incident and scattered fields. The full wave propagation and scattering problem is described by the following two wave equations: where j is a model for the current density on the antenna. Using (9) and subtracting (12) from (11) results in the following expression for the scattered field: We then utilize the Green's function solution of the wave equation to arrive at the Lippmann-Schwinger integral equation [3]: It is important to observe that E sc appears on both sides of the equation. Also note that this equation is nonlinear in that the two unknown quantities, E sc and T , are multiplied on the right-hand side of (14). In order to linearize, and ultimately reconstruct T , we invoke the Born, or single-scattering, approximation. This amounts to replacing E tot with E in in (14). For more details on the Born approximation see [3]. We now obtain the following approximate expression for the scattered field: If we take the Fourier transform of (15) we obtain the frequency domain expression for the scattered field Recall we assume that E in satisfies the scalar wave equation given in (12). In the frequency domain this equation is written as where J is the Fourier transform of j, the current density on the antenna. If we again use the Green's function solution of the wave equation we have where x − x 0 indicates we have taken the unit vector in the direction of x − x 0 . Also note we have used the far-field approximation and assumed that the antenna center is located at the position x 0 . F is known as the radiation pattern of the transmitting antenna and is analogous to the radiation vector [3]. This quantity is proportional to the Fourier transform of the current density. Using this model for the incident field we express the scattered field as Note we have dropped the subscript B for ease of writing. To obtain the expression for the data received we evaluate E sc at the received antenna location, which is x 0 in the mono-static case. Therefore we have where A, the amplitude, is given by where again F t is the radiation pattern of the antenna, and F r is the reception pattern of the antenna. In order to include the antenna motion in our model we assume the antenna follows a path γ ∈ R 3 . This may also be thought of as the flight path for the aircraft which the antenna is mounted on in a SAR system. If we consider a pulsed system where pulses are transmitted at times t n , we have that the antenna has a position γ(t n ). To simplify our analysis we choose to instead parametrize the flight path with a continuous parameter denoted s. We call s the slow time and in contrast now call t fast-time. We utilize these two different time scales because the scale on which the antenna moves is significantly slower than the scale on which electromagnetic waves travel. We therefore have that the position of the antenna at a given slow time s is γ(s). We now replace x 0 with γ(s) in (20) to obtain the following expression for the received mono-static SAR data In the time domain we have Note that for simplicity we will now assume flat topography where x = (x 1 , x 2 ), x = (x, 0). Also observe we have obtained a linear relationship between the data, d, and the target reflectivity T , in terms of the forward operator F. Under mild conditions, the operator F is a Fourier Integral Operator (FIO) [13,4,11,12]. It is well known [4,11,12] that an FIO may be approximately inverted by another FIO. This inversion is called backprojection which we discuss next.

Filtered backprojection.
To form an image we aim to invert (23) by applying an imaging operator K to the collected data. Because F is a Fourier integral operator we can compute an approximate inverse by means of another FIO. This FIO typically has the form of a filtered-backprojection (FBP) operator. In our case, we form an approximate inverse of F as a FBP operator, which first filters the data and then backprojects to obtain the image. We note that for the analysis in this paper we assume we are able to integrate over all of R 2 (i.e. the ground plane) in (22) and again over the full R 2 plane (i.e. the data collection manifold) when forming the image. This is clearly not realistic to any SAR imaging scenario but simplifies the analysis. We will comment on the case of a realistic scenario (i.e. limited ground plane and limited data collection manifold) in the sections that follow. Our imaging operator takes on the form where z = (z 1 , z 2 ), z = (z, 0), D(s, ω) is the Fourier transform of the data in fasttime, and Q is a filter which is determined in several different ways depending on the application. If we insert the equation for the data into above we have In synthetic aperture imaging we are especially interested in identifying singularities, or edges and boundaries of objects, from the scene of interest in our image. The study of singularities is part of microlocal analysis, which also encompasses FIO theory [6,2,12,4]. We define the singular structure of a function by its wavefront set which is the collection of singular points and their associated directions. We use the concepts of microlocal analysis to analyze how singularities in the scene, or in T , correspond to singularities in the image. We can rewrite (26) as where B = KF is known as the image-fidelity operator. Note B, the kernel of B, is called the point-spread function (PSF). We will show that B is a pseudodifferential operator which has a kernel of the form where p must satisfy certain symbol estimates as in the definition for a general FIO [6,13]. It is essential that our image-fidelity operator be a pseudodifferential operator because this class of FIOs has the property that they map wavefront sets in a desirable way. This property is known as the pseudolocal property which states W F (Pf ) ⊆ W F (f ), that is, the operator P does not increase the wavefront set.
In the imaging setting this says that the visible singularities of the scene, T , are put in the correct location with the correct orientation. However, we note that it is possible some edges present in T may not appear in the image, especially if the viewing aperture is limited. Also we observe that singularities in the data due to multiple scattering effects which are not captured by the model may give rise to artifacts in the image.
We now show that the filtered-backprojection operator is indeed a pseudodifferential operator and therefore produces an image where edges are in the correct location and have the correct orientation.
In order to demonstrate that B is the kernel of a pseudodifferential operator our goal is to determine K, the imaging operator, so that B is of the form of P in equation (28). First we must ensure B is an FIO by requiring Q to satisfy a symbol estimate, i.e. we assume for some m Q where K is any compact subset of R × R 2 and for every multi-index α, β, γ, there is a constant C 0 = C 0 (K, α, β, γ).
We now must show that the phase of B can be written in the form Φ(x, z, ξ) = i(x − z) · ξ. In order to determine how close the phase is to that of a pseudodifferential operator, one applies the method of stationary phase to the s and ω integrals. The stationary phase method gives us an approximate formula for the large-parameter behavior of this oscillatory integral. First we introduce a large parameter β by the change of variables ω = βω . The stationary phase theorem tells us that the main contribution to the integral comes from the critical points of the phase, i.e. points satisfying the critical conditions One of the solutions of the above equations is the critical point x = z. Some critical points lead to artifacts in the image (for example the standard right-left ambiguity issue which is well known in SAR imaging [3]), however their presence depends on the measurement geometry and the antenna beam pattern.
In the neighborhood of the point x = z, we use a Taylor expansion of the exponent to force the phase to look like that of a pseudodifferential operator. We make use of the formula where in our case f (z) = 2kr s,z . We then make the Stolt change of variables This change of variables allows us to obtain a new form of the point-spread function where η is the Jacobian resulting from the change of variables, commonly called the Beylkin determinant [2]. From this expression, we see that the phase of B is of the form required of pseudodifferential operators. With the symbol estimate requirements made on Q, we have shown that our image-fidelity operator is approximately a pseudodifferential operator. We conclude that the pseudolocal property holds and therefore our microlocal-analysis-based reconstruction method preserves the singularities or edges in our scene of interest. In order to determine Q precisely we use a standard result in the theory of pseudodifferential operators [11,4]. This result tells us that a pseudodifferential operator can be written as where p(z, ξ) = e −iz·ξ B(e iz·ξ ). The symbol p has an asymptotic expansion where α is a multi-index. In other words, the leading-order term of p(z, ξ) is simply B(z, z, ξ). For more information on this asymptotic expansion of the symbol see [11]. Therefore, according to the symbol calculus, the leading-order term of (27) is given by We now choose the filter Q such that the image is as close as possible to ideal. The ideal image would be of the form Note when we assume the data collection manifold (i.e. the area in ξ space where we have data) is the entire R 2 plane and that we are able to integrate over the full x−plane the above expression equals T (z). It is now clear how we should choose the filter Q given our definition of the ideal image, that is, we choose 4. BLUE calcuation. We now calculate the best linear unbiased estimate of the reflectivity function from the collected data which we assume has the form (22) plus additive measurement noise. Therefore we have that (39) D(s, ω) = e 2ikrs,x A(ω, s, x)T (x)dx + n(s, ω) where we assume n is zero-mean independently, identically distributed noise in s and ω, i.e. we assume where σ 2 is the variance of the noise for a single value of s and ω and δ is the Dirac delta function. We aim to estimate T (x) from measurements D(s, ω) via a linear estimator, therefore we assume the estimate T (z) is of the form We seek to find the weights or filter Q such that the variance is minimized subject to the constraint that E[ T (z)] = T (z) (i.e. the estimator is unbiased). Note this Q is to be determined, it is not assumed to be the filter in equation (38). This amounts to a constrained optimization problem of the following functional, where the first term is the variance of the estimator and the second term includes the unbiased constraint. We will be using the Lagrange multiplier approach from calculus of variations in order to minimize the above functional. We note a similar minimization problem was considered in [15] for the case of imaging a target embedded in pseudostationary clutter.
Lets begin by calculating the variance, i.e. the first term in (42). First note the expected value of the estimator T (z) is where we have used that fact that T, A, Q, and the exponential are deterministic and that E[n] = 0. Now calculating the variance we have We now perform the Stolt change of variables as in the backprojection filter derivation from (s, ω) → ξ, therefore the variance simplifies to where again η is the Jacobian from the change of variables.
Next we consider the second term in the function J (Q) coming from the constraint that the estimator be unbiased. That is, we seek an estimator T (z) such that This is equivalent to the following expression where we have used the definition of T (z) and that fact that E[n] = 0. Given our assumptions on the ground plane and data collection manifold being unlimited we note that this condition is equivalent to where we have used the definition of the Dirac delta function on the right hand side. We note that this is one solution technique to (47) which allows us to make the connection to the backprojection technique. Solving (47) may be thought of as finding Q such that the integral operator's kernel has T as an eigenfunction associated with the eigenvalue one. Other methods of solution of this equation are left as future work. Now we make the Stolt change of variables on the left hand side as we did in the variance calculation. In addition we insert the exponential definition of the Dirac delta function on the right hand side to obtain where we have used the symbol calculus as in the backprojection case to obtain the leading order term of the left hand side. Equation (49) is equivalent to the condition This gives us an expression for the functional J (Q) we wish to minimize To find the minimizer we seek Q such that the variational derivative of J (Q) with respect to Q is zero and that the derivative of J (Q) with respect to λ is also zero. Taking the variational derivative of the functional we obtain d d (J (Q + Q )| =0 = σ 2 Q (z, ξ)Q(z, ξ)η(z, ξ)dξ + σ 2 Q (z, ξ)Q(z, ξ)η(z, ξ)dξ + λ Q (z, ξ)A(z, ξ)e ix·ξ dξ = 2 Re σ 2 Q (z, ξ)Q(z, ξ)η(z, ξ)dξ This expression must hold for all Q . Taking the derivative of J (Q) with respect to λ we obtain dJ dλ = [Q(z, ξ)e ix·ξ A(z, ξ)η(z, ξ) − e i(x−z)·ξ ]dξ = 0 (53) which implies Q must satisfy Substituting this expression for Q back into the variational derivative equation (55) gives us an implicit definition for λ Note this expression for λ need not be solved to determine the BLUE of the target reflectivity, we include it only for completeness. If we insert the expression found for Q into the definition for T (z) we obtain Looking at the first term above we see that we obtain precisely the backprojected image from the previous section, i.e.

5.
Conclusion. In this paper we have investigated the relationship of the backprojection method of reconstructing a target scene from SAR measurements and the best linear unbiased estimate of the target scene. This relationship was studied assuming that we are able to image over the entire ground plane (which is assumed flat), and that the data collection manifold is unlimited (i.e. infinite bandwidth and synthetic aperture). In addition we consider the case of zero-mean white (i.i.d) additive noise. We have found that under these assumptions the best linear unbiased estimate is precisely the backprojected image of the scene when the filter is chosen based on a microlocal criterion. This improves upon the author's previous work in comparing the GLRT and backprojection [14]. We note assuming colored or correlated noise would modify the BLUE solution, in particular the equation for λ (55), however it should not affect the unbiased constraint (49).
We note that in a realistic SAR imaging scenario we are unable to integrate over the full ground plane and the data collection manifold will be limited. In this case the unbiased constraint (47) will not be met given that we utilize the filter determined above. Therefore it is impossible to obtain a truly unbiased estimator by using this backprojection/BLUE estimate in the case of a true SAR imaging scenario. We note that in [16] it is argued that using the backprojection estimate results in an estimate whose first order statistics show the visible singularities (i.e. edges) of the target at the correct location and orientation. They also state that the bias of the backprojection estimator is one degree smoother than the first order statistics of the estimate itself. These results hinge on the fact that the backprojection operator is a pseudodifferential operator. Investigating different solutions to (47) could potentially lead to an estimator that is capable of being truly unbiased in the case of a limited data collection manifold and ground plane. In addition another question to be considered is whether these other solutions would lead to a pseudodifferential operator and hence an estimate with singularities in the correct location with the correct orientation. We leave these questions as future work.