SAR CORRELATION IMAGING AND ANISOTROPIC SCATTERING

. In this paper we investigate the ability of correlation synthetic-aperture radar (SAR) imaging to reconstruct isotropic and anisotropic scatter- ers. SAR correlation imaging was suggested by the author previously in [34]. Correlation imaging algorithms produce an image of a second-order quantity describing an object an interest, for example, the reﬂectivity function squared. In the previous work [34] it was argued that the eﬀects of volume scattering clutter on the image can be minimized by choosing which pairs of collected data to correlate prior to applying a backprojection-type reconstruction algo- rithm. This choice of pairs for the correlation process is determined by what is known as the memory eﬀect of scattering of waves by random scatterers [42, 43, 40, 41, 7, 14]. It is the goal of this current work to determine the dif-ferent imaging outcomes for an isotropic or point scatterer versus an anisotropic or dipole scatterer. In addition we aim to determine if removing contributions to the image due to the memory eﬀect is necessary for diminishing the contri- butions of anisotropic or clutter scatterers to the scene of interest. Finally we extend the analysis of [34] to the polarimetric SAR case to determine whether the additional data provided by this modality contributes to decreasing the eﬀects of clutter on the SAR image.


1.
Introduction. The synthetic-aperture radar modality illuminates a scene of interest with electromagnetic waves and reconstructs an image of this scene from the measurement of the scattered waves using airborne antenna(s). In this work we are interested in understanding the properties of SAR images reconstructed using a correlation imaging algorithm. This algorithm requires first a correlation, or multiplication, of the standard measured SAR data at different values of the slow and fast time parameters. Then a backprojection type reconstruction is performed using this correlated data.
Correlation imaging is of interest because it has been shown to improve resolution, reduce sidelobes, and minimize the effects of clutter scatterers on the final image [34,42,43,40,41,7,14]. Note clutter is a term often used in the SAR literature to describe scatterers present in the imaging scene which are not of interest. It usually refers to non-target scatterers in the scene, such as ground or ocean surfaces and vegetation. Targets, instead, are scatterers in the scene that one is attempting to detect and/or classify such as vehicles and other man-made objects.
Scattering from clutter often obscures the data or scattering information from targets. Therefore it is our goal to remove the effects of clutter (or the presence of clutter altogether) from images in order to better detect and classify targets. Note that clutter in this sense does not describe the propagation medium. In many SAR applications (including the ones considered in this paper) the propagation medium is air which is well modeled by free space. Therefore we are not considering the effects of propagating through a random medium. We also note that all scatterers in this work are assumed to be stationary, i.e. not moving.
The most relevant literature to this current study [42,43,40,41,7,14] seeks to make use of the memory effect of scattering of waves by random media or scatterers to minimize the effects of clutter in electromagnetic imaging applications. The memory effect [16,20] states that the angular correlation function (ACF) of optical waves scattered by random media is only nonzero along a memory line (for rough surface scattering) or memory dots (for volume scattering). These lines/dots correspond to correlations of specific angular pairs of incident and scattered directions.
It is argued in [42,43,40,41,7,14] that this effect holds for electromagnetic waves as well and that by avoiding the inclusion of the memory line or dots in the data one is able to decrease the effect of clutter on an image formed from electromagnetic scattered field data. The angular correlation function is given as follows: wherek ij andk sj are two different incident and scattered field directions for j = 1, 2 and where the hat indicates we are considering a unit vector. F is the scattering amplitude from a given incident direction into a given scattered direction. The angle brackets indicate some sort of averaging is performed. Ideally this averaging is over realizations of the random media, though in most radar imaging scenarios only one data set is available. The memory directions are defined to be the incident and scattered field directions that satisfŷ where we now include the additional subscript x to indicate we are looking at the x component of the incident and scattered field vectors in (1). The subscript indicates we are considering the case when the incident wave is a plane wave with e ikixx dependence, i.e. it is propagating in the x direction. It was shown in [16] that the correlation function in (1) is negligble when you remove the memory directions defined in (2) prior to the calculation.
In [42,43,40,41,7,14] it is argued this method reduces clutter effects on radar images specifically for both volume scatterers and rough surfaces. The author's previous work [34] aimed to study the ability of a correlation imaging algorithm which avoided the memory directions to minimize the effects of vegetation-like clutter in SAR imaging. In this work we aim to expand upon the results of [34] in order to more fully understand the effects of the correlation imaging algorithm on images of isotropic vs. anisotropic scatterers. In particular we will focus on the difference between isotropic and anisotropic scatterers modeled by point scatterers and dipoles respectively. This work expands upon that of [34] to incorporate the case when a dipole is long in comparison to the wavelength. Also we aim to determine whether removing data that corresponds to the memory effect prior to correlation affects the ability of this algorithm to reduce the effects of clutter scatterers well-modeled as anisotropic scatterers or dipoles. In addition we will investigate whether utilizing the polarimetric SAR data modality, which provides multiple data sets and measures the polarization state of electromagnetic waves, provides additional benefits when seeking to minimize anisotropic clutter.
The paper is organized as follows. The behavior of the point scatterer or isotropic scatterer from [34] is reviewed in section 2. In section 3 we discuss the anisotropic scattering behavior in the context of two different scaling scenarios dependent on the length of a dipole scatterer and the resulting correlation image in each scenario. In both sections 2 and 3 we discuss the utility of removing memory directions for each type of scatterer considered. In section 4 we provide numerical simulations to demonstrate the differences between the isotropic and anisotropic analytic image expressions found in section 3. In section 5 and 6 we expand the work in [34] to the polarimetric SAR modality to assess whether adding in polarimetric diversity further lessens the effect of anisotropic scatterers in a correlation backprojection algorithm. In section 7 we compare standard SAR correlation images with polarimetric SAR correlation images for the case of point scatterers embedded in dipole scatterer clutter. We conclude with observations and discussions of future work in section 8.
2. Isotropic scattering. We assume that we observe an area of interest with a monostatic synthetic-aperture radar (SAR). That is, we assume the transmitting and receiving antenna(s) are colocated. We also assume the SAR system measures only one polarization state and therefore we may use a data model resulting from the scalar wave equation. That is, assuming we are in free space, we can simplify Maxwell's equations to the wave equation [9,36] where E is the electric field and c 0 is the speed of light in free space. The full scattering problem can be described by the following 2 wave equations: where E tot = E in + E sc , E in is the transmitted or incident electric field, E sc is the scattered field, c −2 (x) = c −2 0 − T (x), T is called the reflectivity function which models the scatterers on the ground, and j(x, t) is the current density induced on the antenna. The area where scattering takes place is assumed to be a flat surface, i.e. a plane in R 2 .
We assume we are able to measure the scattered field only at the antenna position. The position of the antenna(s) is given by γ(s) ∈ R 3 where γ(s) is the parametrized curve describing the flight path of the antenna platform. This curve is parametrized by s which is known as slow time, which is distinguishable from t, fast-time, as we assume the speed the electromagnetic waves travel is significantly faster than the speed of the platform. We define T (x) to be the reflectivity function or scattering strength describing an object of interest lying at position x on the surface that is being observed. The goal of SAR imaging is to reconstruct, or estimate, T (x) from measured data. Subtracting the second equation in (4) from the first and solving for E sc yields the forward model for the collected data: where A is a (complex) amplitude that includes the antenna beam pattern, the transmitted waveform, geometrical spreading factors, etc. Here ω is angular frequency and R x,s = |R x,s | where R x,s = γ(s) − x is the vector from the antenna to the point on the ground x. F is known as the forward operator. Note this model, which represents the scattered electric field measured at the antenna, is derived in detail in [9] and uses the Born approximation and the start-stop approximation.
In this section we assume we have ideal measurements without clutter contributions. We also assume the effect of matched filtering or correlation reception is applied in the amplitude A and hence system or measurement noise has been filtered out prior to our processing here. In order to reconstruct T (x) from the measurements d(s, t) we seek an approximate inverse of the forward operator F. The inversion algorithm is known as the method of backprojection [9].
2.1. Correlated data. In [34] we developed a correlation SAR imaging algorithm which required a 'correlation' or preprocessing step prior to imaging (or backprojection). Note we put correlation in quotes because we are not taking a statistical average and are simply multiplying two data sets. We have that the correlated data is given by where the overline indicates we are taking the complex conjugate. In the frequency domain we havẽ where we write the correlated data in terms of a forward operatorF and the unknown describing the target C T (x, x ) = T (x)T (x ).
2.2. Correlation imaging. As in [34] we utilize a filtered-backprojection-type algorithm to approximately obtain a reconstruction or estimate of the target function C T . We have from [34] that the image is of the form: where Q is a filter to be determined based on an optimality criterion derived from microlocal analysis. We insert the data expression (7) into the image expression above to obtain where K is known as the the point-spread function (PSF) [9]. Note K is given by (10) and ideally it would be of the form δ(x − z)δ(x − z ). That is, An ideal PSF implies we can reconstruct C T exactly if we are able to integrate over all frequencies and slow time values. This is not practically possible as this requires an infinite amount of data. Therefore it is the goal of filtered backprojection to obtain a point-spread function (PSF) with the phase given above in (11). In this special case the image-fidelity operator (i.e. the composition of the forward and imaging operators) is a pseudodifferential operator, a special class of Fourier integral operator which ensures the edges present in the scene of interest appear at the correct location with the correct orientation in the reconstructed image [36,9,31,22]. This is one method of obtaining an approximately optimal inverse or imaging operator. To that end we begin by applying the stationary phase approximation as in [34] and then performing the Stolt change of variables we obtain where η, η are the Jacobians resulting from the Stolt change of variables. Note this change of variables is given by where Ξ(x, z, s, k) = 1 0 ∇g 1 | x+µ(x−z) dµ and Ξ (x , z , s, k) = 1 0 ∇g 2 | x +µ(x −z ) dµ and g 1 (x) = 2kR x,s and g 2 (x) = 2k R x,s . Using symbol calculus as in [34] we are able to show that the optimal backprojection filter is of the form where χ Ω is a smooth cut-off function that prevents division by zero.

2.3.
Analysis of the point-spread function. With the filter choice above in (14) the point-spread function is approximately of the form An approximate closed form of K, as well as the resolution of the image, is determined by the data collection manifold, which is the region in ξ and ξ − space where we have data [9]. Note that near z = x and z = x we have Ξ(x, z, s, k) ≈ 2kP[R z,s ] and Ξ (x , z , s , k ) ≈ 2k P[R z ,s ] respectively, where P is the operator that projects a vector onto its first two components. We can use this approximation for ξ and ξ to determine an approximation of the data collection manifold. Under certain assumptions, i.e. for a straight flight path and flat topography, the data collection manifold for standard SAR (i.e. ξ alone) consists of a sector of an annulus [9,30]. For narrow angular apertures, that is, narrow antenna beam patterns, the annular sector is closely approximated by a rectangle [9,30]. For the correlation case the data collection manifold, under these assumptions, is a four dimensional parallelepiped whose cross sections are rectangles in the ξ and ξ spaces respectively. We therefore have that (15) can be (approximately) factored: Following the analysis in [34] and assuming a dipole antenna of length L we obtain where ∆k = k max −k min is the bandwidth of the radar system, k 0 = (k min +k max )/2, and sin ψ = z 1 /R where R = z 2 1 + H 2 is the slant range given an altitude H of the antenna platform and sin ψ = z 1 /R. Note ψ is the elevation angle of the flight path when the antenna has the same z 1 coordinate as the point being considered on the ground. We note that in standard SAR imaging we obtain the following PSF It is important to observe that at the peak of the PSF, i.e. when x = z and x = z , the magnitude of the correlation PSF is that of the standard PSF function squared. Hence correlation imaging will increase the magnitude of isotropic scattering strengths in images. In addition, as was discussed in [34], the correlation PSF in the case when z = z has a narrower peak than the standard PSF and reduced sidelobes, furthering indicating the correlation algorithm is preferable for isotropic scatterers.
We now investigate the behavior of K with respect to the memory directions as in [42,43,40,41,33]. In our notation the memory directions are defined by the directions k R z,s and k R z ,s that satisfy: which we denote more compactly with ξ d = 0. To that end we will perform the change of variables given by ξ = ξ + ξ d , where ξ d is equivalent to the definition in equation (19). We will define a simple neighborhood of ξ d = 0 to describe the memory directions we will consider avoiding in the correlation process. We define the neighborhood to be Ω d = ] where 1 , 2 are parameters that may be chosen to eliminate more or less of the incident/scattered field direction pairs. With this change of variables we find the PSF in the neighborhood of the memory directions Ω d is given by We refer to this as the inner PSF. If we choose to remove the memory directions the PSF used in the backprojection algorithm will be given by Note at the peak of the PSF when z = x and z = x we obtain the magnitude which indicates removing memory directions will decrease the magnitude of isotropic scattering strengths. Whether it becomes less than that of the standard SAR image depends on the choice of the memory directions neighborhood Ω d . We will discuss this choice further after analyzing the anisotropic target case.
3. Anisotropic scattering. In order to consider anisotropic, or directional scattering behavior, we now consider a scatterer made up of dipoles. The data model is derived in [34] and we summarize here. We assume that the incident field radiated from the transmitting antenna may be expressed as where F is the radiation vector of the antenna, this includes the waveform and will later become part of the amplitude function A as in equation (5). We then calculate the current on the receiving clutter dipole, following from [9] we have, where L c is the length of the clutter dipole. Note we integrate over all z that lie along the clutter dipole and use the far-field approximation as in [9]. The current induced on the clutter dipole will cause it to reradiate, giving us the expression for the scattered field. We again use the far-field approximation to obtain where F denotes the Fourier transform and J is the current on the antenna. We insert the current found above in equation (24) and calculate the resulting integral. We obtain To obtain the data at the receiving antenna we evaluate the above expression (26) at the antenna location and also include the antenna reception pattern. Assuming our system is monostatic this location is γ(s) which gives us where we have assumed the antenna reception pattern is identical to its radiation pattern, F . Note we assume there is a clutter dipole at every ground location x in the scene of interest. We then integrate over x and our data expression for clutter becomes In addition we may rewrite this using the amplitude function used in (5) as well as the notation R x,s = |x − γ(s)|. We also now include a scattering strength, or reflectivity, ρ c (x), for each clutter dipole, as each dipole may have a different reflectivity. The data received from clutter is therefore expressed as 3.1. Correlated data and correlation imaging. We now consider correlated data from dipole scatterers as in the isotropic case. We havẽ We will apply the same filtered backprojection technique as was used for the point scatterers above. We have that We rewrite this in terms of the Stolt change of variables from equation (13), use the symbol calculus as in the isotropic scattering case, and apply the optimal filter from (14) to obtain We will also investigate the behavior of I C with respect to the memory directions as before by performing the change of variables given by ξ = ξ + ξ d , where ξ d is given in (19). We then obtain the following expression for I C : 3.2. Analysis of the point-spread function. In [34] the behavior of the pointspread function was analyzed. Note returning to the original variables we have that the point-spread function in this case is given by where we note the argument of the sinc 2 functions depends on L c k = 2πL c /λ. This defines a scaling which affects the resulting image.
3.2.1. Short dipole scaling. The simplest scaling scenario is when Lc λ 1. In this case we have Therefore, in this case, the point spread function will have the same form as the isotropic point-spread function modified by a factor of L 4 c . Thus we have We also calculate the inner PSF in the neighborhood of ξ d = 0, Ω d , as was considered in the isotropic case above and in [34]. We obtain the following expression: Now subtracting the above quantity from the full PSF we obtain the outer PSF, or the PSF used when we remove the memory directions, Prior to removing the memory directions we have that the magnitude of the peak of the PSF for the anisotropic scatter is given by and when we remove the memory directions we have a peak value Inverse Problems and Imaging Volume 12, No. 3 (2018), 697-731 We therefore see that in both cases, i.e. with and without the memory directions, the factor L 4 c appears in the anisotropic case and not the isotropic case, and since L c is small we have that the anisotropic PSF will be significantly smaller than the isotropic PSF regardless of whether we remove memory directions or not. We can also compare to the standard SAR PSF for an anisotropic scatterer. This has the form In the case when L c λ we have that the peak value is approximately We note in the standard SAR case the peak value of the isotropic PSF is larger than the anisotropic peak, however there is a greater difference in the correlation SAR case. Hence in this scaling scenario it is desirable to use correlation SAR if the expected clutter has anisotropic behavior, that is, we find that correlation SAR suppresses anisotropic clutter better than standard SAR backprojection.

3.2.2.
Long dipole scaling. We now consider the case when Lc λ ∼ O(1) or larger. In this scenario the sinc 2 functions behave as such and hence the analysis in [34] holds. We quote the relevant results here. Note we have made a correction to the proportionality constant below in this work.
Theorem 3.1. We assume the data collection manifold is closely approximated by a rectangle, i.e. we assume a linear flight path γ(s) = (0, s, h) where h is a constant altitude and s ∈ [s 1 , s n ] ⊂ R and s 1 , s n indicate the start and end values of the slowtime parameter for the specific flight path. We assume the transmitting/receiving antenna is a dipole of length L and that the bandwidth is given by ∆k = k max − k min with center wavenumber k 0 = (k max + k min )/2. We assume flat topography as well as a narrow angular aperture or antenna beam pattern, that is, we assume the angle θ between the vector in the data collection manifold ζ 1 = ξ(k max , s 1 ) − ξ(k min , s 1 ) and ζ n = ξ(k max , s n ) − ξ(k min , s n ) is small. Lastly we assume the orientation of the dipole, e(x) = [ e 1 (x), e 2 (x)], is such that e i (z) = 0 for i = 1, 2 to ensure the denominators below are never equal to zero. In this case the integrals appearing in the clutter point-spread function in equation (34) may be approximately factored. The clutter point-spread function is thus expressed as: In addition in a neighborhood surrounding ξ d = 0 of size 1 × 2 the inner clutter point-spread function has the following form: Note the triangle function tri(x) is defined as follows Also observe that in this theorem we make the somewhat limiting assumption of a narrow antenna beam pattern. This is the same assumption used in standard resolution analysis [9]. This assumption is necessary in order to factor the ξ and ξ integrals in the proof of the theorem. However in simulations it is found that we may have a beam wide enough to see the whole scene of interest at each slow time location and still obtain the results predicted by the analysis.
In [34] the effect of removing the memory directions is also discussed. Note in this case the outer point-spread function is given by It is argued in [34] via a rough estimate that K outer C may be made zero (or as small as desired) via the choice of 1 and 2 . The argument is as follows: first note that the integral present in the outer clutter PSF is actually over the support of the triangle functions present in the integrand of equation (47). The support is the Our aim is to estimate the value of the integral. To that end we maximize the value of the integrand within the (y 1 , p 1 ) plane. The local extremum of the integrand is at the origin when this region of integration is small, i.e. when L c is small. We may therefore estimate This is precisely what we obtain when we assume sinc 2 ∼ 1 in section (3.2.1) .
We now aim to expand these results for the case when L c gets large and the integral estimate above is less accurate. In this case we may approximate the triangle functions in (47) by tri 2 L c e 1 (z) In this case the full anisotropic PSF is given by Since (50) is a product of two identical integrals, we will consider just one integral factor for now. We have Inverse Problems and Imaging where FT is the Fourier transform, * denotes convolution, and in the first line we use the convolution theorem. Also note in the last line we use the assumption that which is true for the simulation parameters used in the following sections. Now using this calculation in the expression (50) for K C we have Next we perform a similar calculation to determine the behavior of K inner C when L c is large. Again we assume tri 2 Lc e1(z) p 1 ∼ 1. In this case we have After calculating the y 1 and p 1 integrals using various Fourier transform results we obtain The details of this calculation are given in the appendix.
Inverse Problems and Imaging Volume 12, No. 3 (2018), 697-731 The peaks of the anisotropic PSFs in this scaling scenario are given by and Recall the value of the peak of the isotropic scatterer's PSF is given by We note that without removing the memory directions the difference between anisotropic and isotropic scatterer PSF peak value depends on the size of 1/L 2 vs. π 4 L 2 c /( e 2 (z) e 2 (z )), which is dependent on the specific parameters of a given imaging scenario. They can be very close to the same value when the anisotropic scatterer is long compared to the wavelength. Hence in this scenario removing the memory directions becomes essential if the goal is to the remove the effect of anisotropic scatterers from the image. This is because it is possible to choose 1 and 2 such that ∆k 2 sin ψ sin ψ − 4π 1 4π∆k sin(ψ ) L − | 1 2 | and hence making the anisotropic PSF negligible in comparison to the isotropic PSF. Note this is partly due to the assumption that 4π/L > ∆k sin ψ. We may also compare the correlation PSF with the standard PSF for the anisotropic scatterer when L c is large. In this case we have which has value at the peak We may recall that the peak value of the isotropic scatterer's standard PSF is given by 16π∆k sin ψ L .
Inverse Problems and Imaging Volume 12, No. 3 (2018), 697-731 It is clear that without removing the memory directions in the case when L c is large it is preferable to use standard SAR when attempting to remove effects from anisotropic scatterers as the peak value of the PSF is much larger in the correlation case for the anisotropic scatterer. If one considers removing memory directions then the correlation algorithm is preferable depending on the parameters of the imaging scenario. We will now demonstrate these results via numerical simulations.
4. Numerical simulations -Isotropic vs. anisotropic. We now simulate monostatic SAR data for a single point scatterer (i.e. isotropic scatterer) and a single dipole scatterer (i.e. anisotropic scatterer) of varying length. We simulate data with simple numerical quadrature using the original forward data model in equation (5). We then reconstruct images of the reflectivity functions of each scatterer type using standard filtered backprojection as well as correlation-based backprojection where we vary the size of the neighborhood of ξ d = 0 removed. The scene on the ground is assumed to be 12 meters by 12 meters represented by 267 by 267 discretized ground locations. We assume the flight path of the antenna is a straight line given by (300, s, 300) meters where −150 ≤ s ≤ 150 with 350 slow time steps. We assume our transmitted frequency band is .5 − 1.5 GHz with 4000 frequency samples (i.e. well above the needed Nyquist sample rate). We assume a wide antenna beam pattern where we are able to see the entire scene of interest at each flight path location. In these simulations we consider a point scatterer located at the center of the scene of interest as well as a single dipole at the center with lengths approximately equal to .3λ, 3λ, and 30λ where λ ≈ .2998 meters is the center wavelength. We define the resolution cell to be approximately of size λ × λ where we therefore have 81 by 81 pixels.
Using the simulated data we form two different SAR images. The standard SAR image is calculated using filtered backprojection as in [9]. Next we use the correlation backprojection as in (8). For the numerical calculation of the correlation image we consider only the case when z = z so the images are two-dimensional as in the standard backprojection case. We begin by rewriting (8) in the time domain. We have, whered is the inverse Fourier transform of the filtered data QD. We first perform the ω and ω integration to obtain where we have factored the data into distinct parts that depend on (s, t) and (s , t ). Next we define t = t +t and s = s +s which yields We recall that in the analysis we used the variables ξ and ξ . These variables arose from the Stolt change of variables from (k, s) → ξ and (k , s ) → ξ and their difference ξ d . In practice we work in the time domain so instead we use (t, s) and (t , s ) and the difference between these twot ands. With that in mind we now perform thet integration to obtain Lastly we integrate over t resulting in the expression we evaluate numerically: Therefore in the numerical simulations we define the memory directions bys = 0. We will exclude these directions by implementing a tolerance on the size ofs, that is we require |s| > tolerance.   We first display the results of a single point target centered in the scene of interest. In Figure (1) we display the results of the standard backprojection vs. the correlation backprojection without removing the memory directions. Clearly the correlation backprojection provides a better resolved image with smaller sidelobes. In Figure (2) we display the results of the correlation algorithm after removing memory directions. We see some slight changes in the resulting image but as demonstrated in Table (4.1) the peak value where the target is located does not degrade much when removing the varying neighborhoods of the memory directions. This is a desired result as we hope the target will be impacted less than clutter, or anisotropic scatterers, when correlation backprojection is used.  1.1491e-04 -6.2788e-06i Correlation BP with memory directions 0.6668e-08 -0.0410e-08i Correlation BP withs = 1 0.6575e-08 -0.0410e-08i Correlation BP withs = 25 0.5030e-08 -0.0377e-08i Correlation BP withs = 50 0.3642e-08 -0.0290e-08i Correlation BP withs = 75 0.2490e-08 -0.0208e-08i Table 2. Peak Image Value at center pixel, single dipole target, orientation π/4, length= .3λ Next we display the results of a single dipole centered in the scene of interest with orientation π/4 and length .3λ. In Figure (3) we display the results of the standard backprojection vs. the correlation backprojection without removing the memory directions. We see the short dipole behaves similarly to the point scatterer and that correlation imaging reduces sidelobes. We also note the much smaller amplitudes due to the L 2 c and L 4 c factors in the standard and correlation images respectively. In Figure (4) we display the results of the correlation algorithm after removing memory directions. We again see some slight changes in the resulting image. In Table (4.2) we see removing memory directions provides a small improvement (i.e. a decrease in peak size) over including them, however it does not have a large impact as expected.   Correlation BP withs = 1 0.6983e-06 -0.1448e-06i Correlation BP withs = 25 0.3869e-06 -0.1183e-06i Correlation BP withs = 50 0.2159e-06 + 0.0199e-06i Correlation BP withs = 75 0.1153e-06 + 0.0139e-06i Table 3. Peak Image Value at center pixel, single dipole target, orientation π/4, length= 3λ We then increase the length of the dipole to 3λ and 30λ in Figures (5), (6), (8), and (9). In these cases we see the anisotropic scatterer behaving as such. We include a second simulated flight path to demonstrate how different the images are Standard BP 0.0090 + 0.0004i Correlation BP with memory directions 0.4206e-04 + 0.0014e-04i Correlation BP withs = 1 0.4072e-04 + 0.0014e-04i Correlation BP withs = 25 0.2535e-04 + 0.0049e-04i Correlation BP withs = 50 0.1363e-04 -0.0087e-04i Correlation BP withs = 75 0.0592e-04 -0.0021e-04i Table 4. Peak Image Value at center pixel, single dipole target, orientation π/4, length= 3λ, second flight path   (7) and (10). The second flight path is given by (30, s, 300) with −50 ≤ s ≤ 50 (all in meters) and N s =150. In the original flight path case the anisotropic scatterer's Standard BP 0.0010 -0.0023i Correlation BP with memory directions 0.3106e-05 -0.0308e-05i Correlation BP withs = 1 0.3053e-05 -0.0308e-05i Correlation BP withs = 25 0.2214e-05 -0.0087e-05i Correlation BP withs = 50 0.1554e-05 -0.0159e-05i Correlation BP withs = 75 0.1115e-05 -0.0061e-05i  Table 6. Peak Image Value at center pixel, single dipole target, orientation π/4, length= 30λ, second flight path image is essentially negligible compared to a standard point scatterer regardless of whether we use standard or correlation backprojection. In the case when the object is 30λ in length the peaks appear at the ends of the object. In the second flight path simulation the radar look direction allows the anisotropic target to be seen fully. In the case when the object is 30λ in length we see that removing memory directions becomes more important and we need to remove a much larger neighborhood of s = 0 in order to mitigate the returns from these scatterers. We can conclude that in the case the anisotropic scatterers are small in comparison to the wavelength that correlation imaging without removing memory directions does a superior job of mitigating anisotropic scatterers in comparison with standard backprojection. As the object increases in length removing memory directions becomes necessary to mitigate these scatterers in certain flight path scenarios. We note this analysis does not take into account the effects of any randomness present in a distribution of clutter scatterers and hence removing memory directions may play a larger role when this additional effect is included in the model. This analysis will appear in a future work.
We also note that all comparison analysis thus far for the isotropic and anisotropic scatterers assumes there is only one scatterer in the scene of interest and therefore comparing peak values of PSFs is applicable. In a true imaging scenario there would be a distribution of clutter scatterers and hence the image would be a sum of many PSFs. We note that every pixel of an image of a distribution of clutter, or anisotropic scatterers, will satisfy the following bound where N is the number of pixels and we assume each pixel is of size λ × λ. Also note max K(z, z, z, z) represents the maximum or peak value of the PSF where the scatterer of interest is actually located. We note that when using the numerical experiment values above this maximal correlation image pixel value (with or without memory directions depending on the dipole's length) is still small enough to ensure we are able to distinguish between clutter (anistropic scatterers) and target (isotropic scatterers) in the case of a distribution of clutter scatterers. If we instead consider standard backprojection the maximal pixel value for a distribution of clutter scatterers is on the same order as the point scatterer pixel value. Therefore in the case of multiple clutter scatterers correlation imaging is preferable. This was shown in the numerical simulations in [34]. 5. Polarimetric isotropic forward model and imaging scheme. Our next goal is to assess whether utilizing a polarimetric SAR system, which gives multiple data sets to correlate, will give improved clutter reduction over standard SAR correlation imaging. We begin by developing a forward model for polarimetric SAR data and then derive the correlation backprojection algorithm for polarimetric SAR. This is of interest because we are better able to model anisotropic scatterers when using the full vector description of electromagnetic waves which polarimetric modalities afford us. 5.1. Polarimetric forward model. We now consider the forward model for the collected data from a monostatic polarimetric synthetic-aperture radar system. We assume the SAR system is made up of two dipole antennas, a and b, which travel along the flight path γ(s) ∈ R 3 . The antennas are assumed to have length L. Note we have that dipole a transmits the waveform p a (t), and the scattered field is received on both a and b. Similarly, dipole b transmits the waveform p b (t), and the scattered field is received on both a and b. The Fourier transforms of the waveforms are denoted by P a and P b . We also assume the dipole antennas have orientations e a = [cos(θ a ), sin(θ a ), 0] and e b = [cos(θ b ), sin(θ b ), 0], respectively. Typically these directions are called H and V, referring to horizontal and vertical linear polarization states.
We model scatterers as a collection of dipoles located at various pixels and with various orientations. We denote by e T (x) = [cos(θ T (x)), sin(θ T (x)), 0] the orientation of the target dipole at location x. Note we will suppress the dependence on x of the target orientation for ease of writing. We assume each target dipole has length L T . We will make assumptions to distinguish between target and clutter scatterers shortly. We aim to develop a forward model in the frequency domain similar to that of the standard SAR forward model: where i = a, b and j = a, b. We say D i,j is the set of the data collected when we transmit on the ith antenna and receive on the jth antenna. Here T is the vector function that describes the target. The "forward" operator F maps the target vector field to the corresponding received data. This operator is derived in [36,35], directly from solutions of Maxwell's equations. It is found in [36,35] that the data when dipole a transmits and dipole b receives is of the form  (a, a), (a, b), (b, a), and (b, b)) as for i = a, b and j = a, b and we define R a x,s × R a x,s × e a := [x a , y a , z a ] and similarly for the b antenna. We denote the amplitude matrix with entries A i,j as A. Also note S(θ T ) is known as the scattering vector of the target and is given by which contains elements of the tensor product e T ⊗ e T which appear after rearranging the vector quantities in (68). In theory for a monostatic system, as seen above, the cross-terms (i.e. sin(θ T ) cos(θ T ) and vice versa) in the scattering vector are repetitive and hence only one is needed. In this case we assume S is a 3 × 1 vector and the amplitude matrix, A, is a 3 × 3 matrix. For simplicity we will assume the target's scattering vector may be replaced with the function below, effectively assuming isotropic scattering by the target (i.e. a point scatterer target): This assumption allows us to define our unknown as ρ(x)S(θ T ) which is independent of the wavenumber k and the flight path position s, as is typical of most SAR imaging problems. This assumption also allows us to find the optimal backprojection filter in a manner similar to [36,35]. We therefore have that the collected polarimetric data vector is given by 5.1.1. Data correlation step. We now perform a correlation operation to the collected measurements as in the standard SAR case. We havẽ where the superscript † indicates we are taking the complex conjugate transpose, ⊗ is the standard tensor product, and vec indicates we are vectorizing the matrix argument by stacking the entries in a single column. We aim to reconstruct the unknown function that we now denote T (x, x ) = vec(T (x)T † (x )). Note we write the data in terms of a forward operator applied to the unknown as in the standard SAR case: whereF is the forward operator.

5.2.
Imaging. We seek to develop a filtered-backprojection-type reconstruction method in order to obtain an image of the quantity T (x, x ). Note this will proceed as in the standard SAR case in [34] and summarized in section 2 above. The filtered-backprojection image will be of the form I(z, z ) = e (−2ikRz,s+2ik R z ,s ) Q(z, s, k, z , s , k )D(s, k, s , k )dkdk dsds (77) where Q is a backprojection filter matrix to be determined below. Plugging in the expression for the correlated dataD into (77) we obtain As before we rewrite the image in terms of the point-spread function K: where K is defined as Similarly to the above scenario the goal is to obtain an image as close as possible to the ideal image that can be obtained from the SAR system in use, i.e. K would consist of delta functions just as in the standard SAR case. We perform stationary phase analysis on the oscillatory integral in (80) as well as the Stolt change of variables as in section 2.
Recall the Stolt change of variables is given by where Ξ(x, z, s, k) = 1 0 ∇g 1 | x+µ(x−z) dµ and Ξ (x , z , s, k) = 1 0 ∇g 2 | x +µ(x −z ) dµ and g 1 (x) = 2kR x,s and g 2 (x) = 2k R x,s . This gives us the following expression for the point-spread function where η and η are the Jacobians resulting from each change of variables (81) respectively.
We then recall a result from the theory of pseudodifferential operators [31,22] which states that an operator of the form (85) can be written as where p(z,ξ) = e −iz·ξ L(e iz·ξ ). We note the symbol p has an asymptotic expansion where α is a multi-index. In simple terms this states the leading-order term of p(z,ξ) is L(z,z,ξ). Here 'leading order' is in the microlocal sense [31], meaning that the higher order terms are smoother than the leading-order ones. To read more about the asymptotic expansion of the symbol see [31]. We therefore have that the leading-order term of (79) is given by It is now clear we should choose the filter Q to be where χ Ω (ξ, z, ξ , z ) is a smooth cut-off function that ensures the numerator above is invertible.
6. Anisotropic case. We now consider the polarimetric SAR correlation image of a dipole scatterer at ground location x with reflectivity ρ C (x), orientation e C (x) = [cos(θ C (x)), sin(θ C (x)), 0], and length L c . Similar to section 3 above the data collected by the radar is given by is the clutter function which is made up of its reflectivity ρ C and its scattering vector which is dependent on the orientation angle of the clutter dipole θ C . Note we will now omit the dependence of θ C on x for ease of notation. We then perform the data correlation step as we did previous cases. Our correlated data is expressed as where C(x, x ) = vec(C(x)C † (x )) and again vec vectorizes its matrix argument.
6.1. Imaging. Next we form an image of the anisotropic scatterer from the correlated data in equation (91). We will use backprojection to form the image just as in the isotropic case. Our image is of the form I(z, z ) = e −2i(kRz,s−k R z ,s ) Q(z, s, k, z , s , k )D C (s, k, s , k )dkdk dsds .
Substituting the expression forD C into the image expression we have Note we will define Q as in (89). In order to analyze I C we first make the Stolt change of variables to ξ and ξ from (s, k) and (s , k ) respectively as in equation (13). We have We recall the definition of ξ and ξ : Recall that near z = x and z = x we have Ξ(x, z, s, k) ≈ 2kP[R z,s ] and Ξ (x , z , s , k ) ≈ 2k P[R z ,s ] respectively, where P is the operator that projects a vector onto its first two components. We see that these approximate values for ξ and ξ appear in the arguments of the sinc 2 functions in equation (93) up to the projection operation. We also note that these vectors appear with a dot product with e C (x), which is assumed to have a zero third component (i.e. the clutter dipoles lie flat on the ground). Thus we may replaceR z ,s with P[R z ,s ] and substitute k R x,s with ξ in the arguments of each sinc 2 in (94).
Next using the symbol calculus and the filter definition as in equations (88) and (89) we obtain where I is the identity matrix. We now introduce ξ d for which ξ d = 0 defines the memory directions. That is, we again let ξ = ξ + ξ d , where ξ d is equivalent to the definition in equation (19). Performing the change of variables we obtain the following expression for I C : We now consider the point-spread function as we did in the scalar SAR case. We have that the clutter PSF is given by Note that the point-spread function is precisely the same as in the scalar SAR work in [34] and in previous sections. We therefore obtain the same analysis as in the previous work and sections above.
7. Polarimetric vs. standard SAR correlation imaging -Numerical simulations. We now simulate monostatic polarimetric SAR data and obtain images using our correlation backprojection algorithm and compare with images obtained using standard polarimetric backprojection. The data is simulated simply by implementing (74) and (90) numerically for the target and clutter scatterers respectively. In addition we will calculate and compare the signal-to-clutter ratio for all the various algorithms considered here.
Note we assume our polarimetric SAR system is made up of two linearly polarized dipole antennas, i.e. H and V polarizations, which are colocated. The scene on the ground is assumed to be 24 meters by 24 meters, represented by 81 by 81 pixels. That is, our resolution cell size is .29 meters by .29 meters. The endpoints of the scene are (±12, ±12) (meters). We assume we have two point scatterer targets embedded in dipole clutter. We assume a clutter dipole is present in every pixel. We consider uncorrelated clutter for simplicity. The orientation of each dipole is assumed to be θ C = π/4 and hence e C (x) = [ √ 2/2, √ 2/2, 0] T . The clutter reflectivity values ρ C (x) are modeled as i.i.d complex Gaussian random variables with zero mean and unit variance. We collect simulated data for a linear flight path which has constant altitude of 300m, the x 1 position is fixed at 300m, and the x 2 position varies linearly between −150m and 150m with 300 slow time sample points. The frequency band used for transmission is 1 − 1.5GHz with 4000 frequency steps.
Using the simulated data we form two different SAR images. First we use the correlation backprojection as in (77). For the numerical calculation of the correlation image we consider only the case when z = z so our images for each polarization channel are two-dimensional as in the standard backprojection case. The algorithm will be the same as given in the standard SAR case, i.e. as in equation (65) however the data will be a vector quantity in the polarimetric case. Again we will define the memory directions bys = 0. We will exclude these directions by implementing a tolerance on the size ofs, that is we require |s| > tolerance.
We compare the correlation images to images formed using a standard polarimetric backprojection algorithm. If we return to the forward data model for polarimetric SAR prior to any correlation processing we have, where we utilize a standard backprojection filter given by where χ Ω is a smooth cut-off function that ensures the numerator is invertible and η is the Jacobian from the Stolt change of variables. This is similar to the backprojection algorithm derived in [36,35]. Figure 11. Two point targets, dipole clutter, each dipole with orientation π/4, length= λ, and Gaussian reflectivity; standard polarimetric reconstruction, HH, HV, VV images respectively In Figure (11) we display the results of standard polarimetric SAR backprojection for the two point targets embedded in dipole clutter in the case when data signalto-clutter (SCR) ratio is 0 dB. In addition we display the results of the correlation backprojection in Figures (12) and (13) fors = 1, 50 slow-time steps respectively. Note we display only the correlations of HH − HH, HV − HV , and V V − V V , that is, each polarimetric data set correlated with itself. The algorithm produces 9 images in total but for brevity we display only these three, the remaining images are all similar qualitatively. We see that the correlation algorithm removes many of the effects of the dipole clutter and improves resolution just as in the scalar SAR case shown in [34]. We note increasing the size of the neighborhood around the memory directions makes no significant difference in the image. This is expected because the size of each dipole is on the order of the transmitting wavelength. In addition we did not add any statistical correlation between the various dipoles. The significance of correlation will be addressed in a future work.
Next we display a series of tables with the image signal-to-clutter (SCR) ratios for each algorithm considered in this paper. Note we consider the cases of data signal-to-clutter ratio equal to 20, 10, 0, −10, −20 dB. We implement this level of clutter by adding the data produced by the two point targets to the dipole data multiplied by the following factor: scale factor = s,k (|T HH (s, k)| 2 + |T HV (s, k)| 2 + |T V V (s, k)| 2 ) s,k (|C HH (s, k)| 2 + |C HV (s, k) where T HH is the data collected from the target scatterers by transmitting on the H antenna and receiving on the H antenna, and similarly for T HV , T V V . Also C HH is the data collected from the clutter scatterers by transmitting on the H antenna and receiving on the H antenna, and similarly for C HV , C V V . Note for the case of scalar SAR we would simply have one term inside each summation. We calculate the image SCR via the following formula SCR = 20 * log 10 where the above norm is the L 2 norm and we repeat this calculation for each imaging algorithm and each data channel in the polarimetric cases.
In Table ( (100). Note here we display the results for the various antenna pairings available with a polarimetric SAR system. That is, we include the cases when we transmit on the horizontally polarized dipole and receive on the same one, transmit on H and receive on V, and finally transmit on V and receive on V.
Lastly Tables (7.3)-(7.5) give the results of the polarimetric correlation algorithm for tolerancess = 1, 10, 50 slow-time steps respectively. We find expected results here that mimic the analysis of the polarimetric algorithm. As expected adding polarization diversity does not increase the image SCR when compared with scalar SAR. Of course the correlation algorithms perform significantly better. This is expected as we saw in the direct comparison of isotropic to anisotropic scatterers, the magnitude of the image of isotropic scatterers is much larger than that of the anisotropic dipole scatterers. We conclude again that correlation imaging does minimize the effects of dipole clutter on SAR images, however simply scalar SAR is needed to achieve this benefit. In addition in the scenarios studied here removing the memory directions is not necessary to see large reductions in clutter. Input Table 11. Input vs. ouput SCR in dB, pol SAR correlation BP, two point targets in dipole clutter,s = 50 scatterers and found that correlation imaging increases the resulting image magnitude for a point scatterer and decreases the magnitude for a dipole scatterer when compared to standard SAR backprojection images. Therefore we validate the results of [34] which stated correlation SAR provided improved clutter reduction in comparison to standard SAR when imaging point targets embedded in dipole clutter. We extended the results of [34] analytically for two scaling scenarios, i.e. when the dipole is small in size compared to the transmitted wavelength as well as when it is comparable or larger in size. We found correlation imaging was superior in both scalings. We also considered whether removing memory directions prior to correlation was necessary to reduce the magnitude of the image from anisotropic scatterers. It is found that while there is some reduction when memory directions are removed, it is not significant enough to warrant the additional consideration of the size of the neighborhood of the memory directions to remove prior to imaging. However, the greatest effect was found in the case when the dipole was longer than the transmitted wavelength.
Lastly we considered whether adding polarimetric diversity to the imaging scenario prior to correlation imaging would improve clutter mitigation. It is found that while correlation polarimetric SAR imaging is preferable to standard polarimetric SAR backprojection in terms of clutter mitigation it does not provide significant benefits over using standard scalar SAR. We conclude by noting this work leaves the important consideration of the randomness of clutter to a future work. × sinc 1 (z 1 − x 1 − y 1 ) sinc 2 z 2 − x 2 − e 2 (z ) e 1 (z ) y 1 × sinc ∆k sin(ψ) ) − e 2 (z ) e 1 (z ) y 1 − e 2 (z) e 1 (z) p 1 dy 1 dp 1 .

(107)
Now in order to find the y 1 integral in (105) we must convolve the results of the two calculations above and evaluate the result at 2k 0 sin ψ 1 − e2(z ) e1(z ) e1(z) e2(z) . We arrive at the following Now combining the results of the y 1 and p 1 integrations we obtain that when L c is large K inner C is approximately given by