Microlocal analysis of Doppler Synthetic Aperture Radar

We study the existence and suppression of artifacts for a Doppler-based Synthetic Aperture Radar (DSAR) system. The idealized air- or space-borne system transmits a continuous wave at a fixed frequency and a co-located receiver measures the resulting scattered waves; a windowed Fourier transform then converts the raw data into a function of two variables: slow time and frequency. Under simplifying assumptions, we analyze the linearized forward scattering map and the feasibility of inverting it via filtered backprojection, using techniques of microlocal analysis which robustly describe how sharp features in the target appear in the data. For DSAR with a straight flight path, there is, as with conventional SAR, a left-right ambiguity artifact in the DSAR image, which can be avoided via beam forming to the left or right. For a circular flight path, the artifact has a more complicated structure, but filtering out echoes coming from straight ahead or behind the transceiver, as well as those outside a critical range, allows one to obtain an artifact-free image. Initially derived under a start-stop approximation widely used in range-based SAR, we show that some of these results are robust and hold under a more realistic approximation.


Introduction
Synthetic aperture radar (SAR) systems transmit electromagnetic waves from an airborne or spaceborne antenna, which then scatter from objects of interest on the terrain, e.g., roads, vehicles, buildings or natural features.The scattered waves are measured by a receiver either co-located with the transmitting antenna (monostatic SAR) or on one or more other platforms (bi-or multi-static).Radar transmitters operate in a variety of modes, with the emitted waves ranging from wideband pulses to ultra-narrowband, continuous wave (CW) signals.The Doppler Synthetic Aperture Radar (DSAR) system we analyze in this paper is of the latter type, with narrow temporal windowing superimposed on the received scattered waves generated by a single frequency transmitted wave.(This approach to applying temporal windows to CW signals has some similarity to what is known as pseudo-pulse processing in the radar literature, which has been used previously for SAR [34].)Possible uses of DSAR include low power applications and imaging through media which are either dispersive or have frequency-dependent attenuation.
The narrowband nature of the DSAR waveforms allows accurate measurement of Doppler frequency shifts, which in turn can be used to obtain accurate measurements of relative velocity.Knowledge of relative velocity between stationary scattering objects and the antenna, whose trajectory is assumed known, is then used to locate and image the objects.The dual concept of using Doppler shifts to image rotating or moving objects, using measurements by a stationary antenna, was proposed some time ago by Thomson and Ponsonby [38]; see also [35].More recently, the use of a moving antenna over a horizontal terrain was proposed in Borden and Cheney [5].In Wang and Yazici [40], this approach was named Doppler Synthetic Aperture Radar and examined for bistatic data acquisition; related ideas were further developed in [9,37,42,43,41,44].The reader can also consult these papers for more extensive surveys of the earlier literature.
In contrast to DSAR, in standard SAR the transmitted fields often consist of a train of pulses of short temporal duration and complicated spectrum; the short duration allows for high-resolution estimates of the time between transmitting and receiving antennas (whether in monostatic, bi-static or multi-static geometries), and this travel time is used to estimate the distance (range) between the antenna and the scattering objects.This process is repeated from many locations along the transmitter flight path, providing distance estimates from many different viewing angles, from which both the range and the along-track position of scattering objects can be determined.
In this paper, we study monostatic DSAR using methods of microlocal analysis.Microlocal techniques allow one to develop backprojection algorithms for reconstruction of features on the ground (which we call a scene) up to smooth errors, as well as analyze obstructions to such reconstruction in the form of imaging artifacts (fictitious features in reconstructions of a scene).Such methods have been highly successful in either revealing artifacts, arising in conventional monstatic and bistatic SAR, or in rigorously explaining artifacts that have already been observed.See [12], [30] and [36] for examples of artifacts that are either predicted or explained by microlocal analysis.Also see [13] for an example of how to ameliorate the effects of artifacts.
We describe what types of artifacts appear in DSAR and obtain sufficient conditions under which they can be excluded, allowing for robust imaging of sharp features on the ground (such as walls and edges) via filtered backprojection.
For simplicity, we treat the problem of imaging a stationary scene located on a flat terrain.After considering the structure of DSAR for general flight paths, we focus on two model trajectories, namely a straight flight path and a circular one.For these we characterize what artifacts arise in the imaging, and describe criteria for avoiding them (see Sections 6 and 7).
More precisely, in this paper we first use the Born approximation and other approximations to create a forward map F taking the scene V (x) to be imaged to a windowed Fourier transform W (s, ω) of the DSAR data.(Here, x are coordinates on the ground, s is the slow time parametrizing the transceiver's flight path γ(s) (assumed smooth) and ω denotes frequency.)We then analyze F, showing that it is a Fourier integral operator, and study the implications of its microlocal geometry for imaging by filtered backprojection.
For the linear trajectory, we show in Theorem 2 that there is a left-right artifact, similar to that which arises in monostatic SAR (see Remark 2), and which can be avoided using beam forming to the left or right of the flight track.Furthermore, it is shown in Sec. 8 that this conclusion still holds under a more refined and realistic model, indicating that our conclusions are robust.For a circular flight path, the geometry is different and the artifacts are more complicated (see Remark 7), but in Theorem 3 we are able to precisely characterize the artifacts, and then give criteria for avoiding them (e.g., see (64) in Lemma 5).

The Doppler SAR model
We consider an aircraft-or satellite-borne antenna, transmitting a time-harmonic signal e −iω 0 t that scatters off the terrain and is then detected with the same antenna.The received signal is then used to produce an image of the terrain; we follow [5] closely for the initial modeling; in particular we ignore polarization and use the scalar, time-domain wave equation, (1) where y ∈ R 3 , E is (one component of the) electric field, f describes the source, and the function c is the wave propagation speed.
For the majority of this paper, we take the source f to be of constant angular frequency, of the form f (t, y) = e −iω 0 t δ(y − γ(t)), where γ := {γ(t) | t min < t < t max } is the curve describing the antenna flight path.In particular, in this model we assume that the antenna radiates isotropically.More generally, following [8] and [31], we recall below how an explicit antenna beam pattern can be incorporated in f , but this will only affect the amplitudes and not the phase function and geometry.
To avoid irrelevant degeneracies, we impose some assumptions, Assumption 1.The flight path is strictly above the ground and, in the region between the flight path and the ground, c(y) = c 0 , the constant speed of light in dry air.
We define the reflectivity function, 1 , which encodes changes in the propagation speed.Electromagnetic waves attenuate rapidly as they propagate into the earth, which is modeled by the following assumption: We provide below a very brief sketch of the linearized model of the scattered waves and we refer the reader to [31] for more details.If we formally linearize (1) by writing c = c 0 + δc, E = E 0 + δE, then the scattered field, δE, approximately satisfies (2) where the incident field, E 0 , satisfies (1) with c replaced by c 0 .This is the Born approximation of the scattered field.If we wish to include an antenna beam pattern in (1), we set f (t, y) = J s (y)e −iω 0 t , where J s (y) is related to the current distribution on the antenna.Finally, for the time being, we make the following simplifying assumption: Start-stop approximation: The speed of wave propagation is so large, relative to the motion of the transceiver along γ, that the point where the wave is detected by the receiver after scattering off the terrain is the same as where the transmitter emitted it.
The start-stop approximation is common in range-based SAR.There, the emitted wave consists of a chain of short duration pulses and, at least for airborne systems, the speed of the platform is small relative to the speed of light; see [39,7] for further discussion.In our setup, the start-stop approximation at first glance seems much harder to justify: rather than using short duration pulses (typically with complex spectra), it uses a long duration/narrow bandwidth signal, idealized as continuous wave (CW).We make this assumption because the calculations needed for the microlocal techniques used become far more complicated without it.However, we are able to partially justify the start-stop approximation for DSAR, ex post facto, by showing that the canonical relations associated to the linearized forward scattering map are either structurally stable (in the case of a circular flight path) or robust to a first-order correction (in the case of a linear flight path, as shown in Sec. 8).
For The start-stop approximation yields the total time of travel for a transmitted wave, T tot , as follows: the downward travel time of the wave emitted at time t tr = t from γ(t) to x is c −1 0 R(t); the wave is thus incident to x at time t in = t + c −1 0 R(t).The scattered wave is then detected by the receiver at time t sc = t + T tot , with T tot determined implicitly by ( 4) Under the start-stop approximation, we have (This calculation will be refined in Sec.8.) Using the assumptions above and the Green's function for the wave equation, ( 5) convolving G with f to get E 0 and then convolving G with the right-hand side of (2), we see that under the stop-start approximation, the scattered wave measured on the antenna is The amplitude p(x, ω) is related to J s and is the antenna beam pattern, referring to the fact that one can phase the individual elements of the antenna to interfere in such a way as to selectively illuminate desired parts of the scene on the ground.Note that a further refinement is also possible with regard to beam-forming on the received signal [31].For the majority of our discussion, the beam pattern is not important and unless otherwise stated, we take p ≡ 1 in (6).Note also that, for ease of notation, we suppress the dependence of R and R on x.Furthermore, the slowly varying range factor in the denominator of (6), as well as any constants, will be absorbed into the amplitude a in (10) below.
The Doppler problem, like most radar problems, involves multiple time scales.The speed of light is c 0 ≈ 3 • 10 8 m/sec, whereas aircraft speeds are typically subsonic, or less than about 3 • 10 2 m/sec.Even satellites in low earth orbit travel only on the order of 8 km/sec ≈ 10 4 m/sec.Moreover, the frequencies involved are usually very large, typically above 1 GHz = 10 9 Hertz.To analyze the signal (6), we introduce a "slow time" s, and multiply the data by a windowing function, whose duration is (i) small relative to the antenna motion (i.e., the distance the antenna travels during the window is small relative to a wavelength), but (ii) large enough so that the transmitted signal undergoes a sufficiently large number of cycles over the support of the window so as to be amenable to Fourier analysis.We take the window about t = s to be of the form (ω 0 (t − s)), where (t) is smooth, identically equal to 1 for |t| ≤ L and supported in |t| ≤ 2L, for some appropriately chosen L > 0. Although of compact support, it is natural to consider (•) as being a symbol of order zero (see below), as the functions (ω 0 •) are symbols of order 0 uniformly in ω 0 and L.
Computing the resulting windowed Fourier transform of the data, localizing near s in the time variable and using ω for the frequency variable, we form In (7) we Taylor expand R(t) about t = s (cf.[5, Eqns.(7-8)]) as ( 8) using the convention throughout that dot denotes differentiation with respect to time.Keeping only the linear terms then results in the approximation of W 0 which will be the object of study: After the change of variables t → τ = t − s, this becomes We define the DSAR transform by F : V (x) → W (s, ω), and show below that from the form (10) one can identify F as a Fourier integral operator (FIO).Although not stated explicitly in (10), we will implicitly assume that W has been multiplied by a smooth cut-off function with compact support.This makes physical sense but also is necessary to avoid artifacts in the backprojected image later on.We will describe the main ideas and results on FIOs and other aspects of microlocal analysis as needed, but also refer the reader to [10,23,25] for more detailed accounts.

The DSAR transform as a Fourier Integral Operator
A Fourier integral operator (FIO) is an integral operator whose kernel has the form K(y, x) = e iφ(y,x,τ ) a(y, x; τ )dτ , where the phase φ and amplitude a satisfy certain conditions outlined below.In our case, the output variables are y = (s, ω); the input variables are the components of x; a is as in (10); and φ(y, x) = τ (ω − ω 0 + 2ω 0 Ṙ/c 0 ).Theorem 1.Under assumptions 1 and 2, the mapping F is a Fourier integral operator of order −1/2 associated with the canonical relation C in (14) below.
Proof.To check that F given by ( 10) is an FIO, we need to check the following conditions on the amplitude and phase of (10).
First, the phase φ(s, ω, x; τ ) := τ (ω−ω 0 +2ω 0 Ṙ(s)/c 0 ) is an operator phase function in the sense of Hörmander [23], meaning that its differential in all the variables, (11) is nonzero for τ = 0.This is easily checked.Second, the amplitude a(s, ω, x; τ ) must satisfy certain symbol estimates.We note first that the function (ω 0 τ ) is a symbol of order 0 in the phase variable τ , which means that it belongs to the class x 3 and N denotes the nonnegative integers.Moreover, since the amplitude a is a product of a smooth, nonzero function of (s, x), independent of ω and τ , with the order 0 symbol (ω 0 τ ), a is also a symbol of order 0, i.e., its derivatives satisfy (12) with (ω 0 τ ) replaced by a(s, ω, x; τ ).
Note that, since (•) is of compact support, so is a (as a function of τ ).Thus, F is strictly speaking a smoothing operator; however, in order to understand how singularities in the reflectivity function are transformed into the data, we consider a as a symbol of order 0, which is uniform in ω 0 and L, as mentioned above.Note also that a is of compact support in the spatial variables, since we assume that the antenna is compactly supported.
Finally, the order of the operator comes from the Hörmander convention for the orders of FIOs [23]: This completes the proof of Theorem 1.
The properties of an FIO involve the geometry of certain key sets.The first is called the critical manifold of φ, which for F is defined by Because φ is nondegenerate in the sense that the full gradient of d τ φ in all the variables is nonzero, it follows that Crit φ is a smooth, codimension one surface.
In general, a nondegenerate phase function determines, in addition to Crit φ , a second key set, C, called the (twisted) canonical relation of φ (or F), which is a smooth, immersed submanifold of the total cotangent space of all the variables, both input and output.The canonical relation describes microlocally how F transforms the locations and directions of singularities in the input variables (i.e., of the reflectivity V ) to singularities in the output variables (i.e., of the signal d).Examples of canonical relations include, but are not limited to, the graphs of canonical transformations between the cotangent spaces of the the input and output variables.
The twisted canonical relation C is the image of Crit φ under the map Here T * R 2 denotes the cotangent bundle of R 2 , which for the purpose of calculations can be considered as simply R 4 .With respect to the coordinates (s, ω, σ, Ω; x, ξ) on Here, 0 denotes the zero section of T * R 2 , namely {(σ, Ω) = (0, 0)} and ξ = {(0, 0)}, resp.Since both ω 0 = 0 and τ = 0, the claim above that C does not intersect the zero section of either factor space follows from Remark 1 in Sec. 5 below.From ( 14), one sees that x, s, τ are coordinates on C, which we will use for calculations.
In summary, F is an FIO of order −1/2 associated with C, denoted by In later sections, we study the projections from the canonical relation C to the two factor spaces T * R 2 on the left and right.The left projection Understanding these projections gives us information about the geometry of C and thus the properties of F. The next subsection gives the background needed for the subsequent analysis.
3.1.The left and right projections and the Bolker condition.It is a basic aspect of microlocal analysis that, for a general canonical relation, the singularities of the projections π L and π R , in the sense of C ∞ singularity theory [18], e.g., points where their differentials have less than maximal rank, have important implications for the operator theory of associated FIOs.In particular, in applications to inverse problems via the study of linearized forward maps, the singularities of π L and π R determine whether reconstruction via filtered backprojection (modulo smooth errors) is possible, and allows the characterization of artifacts cf.[32,11,15,17,1,14,33,16,2].For any canonical relation, say C 0 ⊂ T * R n × T * R n associated with an FIO F 0 , if one of the two maps Dπ L and Dπ R is nonsingular at a point λ ∈ C 0 , then so is the other.These derivatives Dπ L and Dπ R are each represented by a (2n) × (2n) matrix, so the nonsingularity condition takes either of the forms, (17) det(Dπ L )(λ) = 0 if and only if det(Dπ R )(λ) = 0.
We denote by Σ the set If Σ = ∅, we say that C 0 is nondegenerate; otherwise, it is degenerate.If the complement of Σ contains a sufficiently small neighborhood of λ 0 , i.e. the determinant of ( 17) is nonzero in a neighborhood of λ 0 , then C 0 is a local canonical graph near λ 0 = (y 0 , η 0 , x 0 , ξ 0 ) in the sense that C 0 is the graph of a canonical transformation χ : T * R n → T * R n defined near (x 0 , ξ 0 ).See, e.g., [23], [24,Thm. 21.2.14].
Even though Dπ L and Dπ R drop rank on the same set, and their kernels have the same dimension, the maps π L and π R may have different types of singularities.
Backprojection methods attempt to form an image by applying the adjoint F * 0 to the data F 0 V .If C 0 is a local canonical graph, the formation of the composition F * 0 F 0 is covered by the transverse intersection calculus for FIOs [23,25,10], resulting in If the canonical relation D contains only points of the diagonal ∆, then F * 0 F 0 is a pseudodifferential operator.In this case, under an illumination assumption, one can construct a left parametrix for F * 0 F 0 , i.e., a pseudodifferential operator Q of order −2m satisfying QF * 0 F 0 = I up to a smoothing operator ( assuming that F * F is elliptic).This means that QF * 0 is a filtered backprojection operator that allows us to reconstruct a function V (x) from the data F 0 V , up to a smooth error.Any sharp features in V (x) (such as discontinuities or edges) that are visible in the data will be present in the image, and vice versa; we say that there are no artifacts in the reconstruction.
If, on the other hand, D contains off-diagonal points, not in the diagonal ∆, using straightforward backprojection reconstruction results in artifacts, i.e., spurious features in the image which are not present in the original scene.This is prevented if the Bolker condition [20] is satisfied, which in the DSAR setting reduces to (i) Σ = ∅ (i.e., the projections are nonsingular everywhere); (ii) π L is one-to-one (injective).(18) An injectivity condition on π L such as (ii) is natural for the prevention of artifacts, since injectivity implies that different points in the scene correspond (microlocally) to different points in the data, while (i) implies that the standard transverse intersection calculus of FIO applies to the composition F * 0 F 0 .When combined into the Bolker condition (18), these ensure that D ⊂ ∆, that F * 0 F 0 is a pseudodifferential operator, and that there are no artifacts in reconstructions made using F * 0 .

Contributions of this paper
We have shown above that the DSAR operator F is an FIO and consequently can be analyzed by the techniques of microlocal analysis.In the remainder of this paper, we analyze the geometry of the canonical relation C, and its effect on the presence (or absence) of artifacts in backprojected images, by studying the geometry of C, and its implications for reconstruction of the scene V (x) via filtered backprojection for two model geometries: when the transceiver trajectory is either a straight line or a circle, with flight path at constant altitude over flat topography.
For the straight-line trajectory, we show in §6 that C is degenerate over the projection of the flight path onto the Earth's surface, where it has a fold/blowdown degeneracy, i.e., π L and π R have singularities of (Whitney) fold and blowdown types resp.; descriptions of these singularity classes can be found in Sec.5.2.Our main result on characterization of F for the straight flight path is Theorem 2. If the flight path γ is a straight, horizontal line, then the DSAR transform F is an FIO of order −1/2, associated to a canonical relation C with a fold/blowdown degeneracy.There is a left-right artifact about the projection, Γ of the flight path γ onto the ground.By suitable beam forming, the artifact can be eliminated, in which case the data FV determines any scene V supported away from Γ, up to a possible C ∞ smooth error.
The proof of Theorem 2 will be given in §6.
Canonical relations with fold/blowdown singularities were studied by Guillemin [21] and, in the context of standard monostatic SAR with a straight flight path and an isotropic antenna pattern, by Nolan and Cheney [31] and Felea [11] 1 .Such canonical relations correspond to problems in which naive filtered backprojection results in left-right artifacts, with objects to one side of the flight path appearing in the image on both sides.Consequently, whenever it is possible, such systems use side-looking antennas, so that only one of the ambiguous locations is illuminated.With a side-looking antenna beam pattern, the Bolker condition is satisfied (which, in this equi-dimensional setting means that C is the graph of a canonical transformation χ : T * R 2 → T * R 2 ), and consequently artifact-free reconstruction of the scene V (x) from the data W (s, ω) is possible.This allows for stable reconstruction of the scene V (x) from the data W (s, ω) by filtered backprojection, without any artifacts due to geometry.
On the other hand, for a circular flight path we show in the lengthier analysis in §7 that the forward operator F has a canonical relation C with a more complicated fold/cusp degeneracy (i.e., π L has a fold singularity and π R has a cusp singularity).There is no longer a simple left-right artifact; it is a singular canonical relation, called an open-umbrella (See Remark 7).However, by an appropriate choice of antenna beam pattern, one can again restrict the microlocal support so that the Bolker condition is satisfied, and this enables stable reconstruction of the scene, up to a smooth error.We give explicit criteria for portions of the scene that can be imaged in this way.Our main result for the circular flight path is the following, the proof of which is given in §7.
Theorem 3. If the flight path is circular, the DSAR map F is an FIO of order −1/2 associated to a canonical relation C with a fold/cusp degeneracy.For a scene V with suitable support, or by suitable beam forming, the associated artifact can be eliminated, and then the data FV determines V up to a possible smooth C ∞ error.
In §8 we modify the start-stop approximation used in the main analysis by adding a first-order correction term, modeling a nonzero transit time for the wave.As a result, the scattered wave is being received at a later time, and thus a different point along the flight path, than when and from where it was transmitted.We show that the result in Thm. 2 characterizing F and its canonical relation for the linear flight path under the start-stop approximation is stable with respect to this correction, even though the blowdown singularity class to which π R belongs is not structurally stable and thus is sensitive to general perturbations.This supports the robustness of the results derived under the start-stop approximation.
We begin, in the next section, with generalities needed for the analysis of the straight and circular flight paths, but which would be applicable to other geometries as well.

Notation, key properties, projections and singularities
We will need properties of the range function R defined in (3) and its derivatives.Recall from (3) that the range vector is R(s) := (x, 0) − γ(s) and the range function is R(s) := |R(s)| (with x suppressed in the notation).Denote the unit vector R(s)/R(s) by R. Let J denote the natural inclusion of R 2 into R 3 , J(x 1 , x 2 ) := (x 1 , x 2 , 0), and J * its transpose, which is the projection J * (ξ 1 , ξ 2 , ξ 3 ) := (ξ 1 , ξ 2 ).Finally, we also define the scaled projection, (19) P With the convention throughout that the dot above a symbol means partial differentiation with respect to s, one easily verifies the following identities: Remark 1. Assumption 1 implies that J * P γ = 0, since R is never horizontal, and therefore from ( 14) and (21) we see that ξ = 0 in C.

5.1.
Properties of the right and left projections.The right projection ( 16) .
Recalling that ξ = −2τ ω 0 d x Ṙ/c 0 , we have (using the parametrization ( 14)) We will thus examine conditions under which d x R and d x Ṙ are linearly independent and, when they are not, analyze the kernel of (22).
This in turn reduces to the question of the injectivity of the map (25) x → ( Ṙ, R).
For two points x, y ∈ R 2 , labeling the associated quantities with subscripts for clarity, the condition Ṙx = Ṙy is easily understood as the Doppler condition.By (20), it says that the down-range relative velocity Ṙx = − R x • γ to x is the same as that to y, or alternatively that the unit vectors R x and R y lie on the same circle in the plane perpendicular to γ, or alternatively, that x and y must lie on the same circular cone with vertex γ and axis γ.On the other hand, the condition Rx = Ry seems harder to characterize.
In the next two sections we focus on two particular flight trajectories, either straight or circular.For simplicity, as noted above we will assume that the flight path is at a constant altitude over a flat landscape, although that is not necessary for the application of the microlocal approach.5.2.Singularities of smooth maps.For the convenience of the reader, we summarize the needed definitions of singularity classes; see, e.g., [4,6,18] for more detailed treatments of singularity theory of smooth functions.Let M and N be manifolds of dimension n; f : N → M be a C ∞ function; and define Σ ⊂ M , For the classes considered, we make the basic assumption that d (det (Df (x))) = 0, so that (i) Σ is a smooth hypersurface, and (ii) at points of Σ, dim (ker Df ) = 1.There then exists a (nonunique) kernel vector field, i.e., a nonzero vector field V along Σ such that Df (p) (V (p)) = 0 for all p ∈ Σ.A map f satisfying these conditions is said to have a corank one singularity.We now define some classes of corank one singularities.
Definition 1. f is said to have a (Whitney) fold singularity along Σ if it only has corank one singularities and, in addition, for every p ∈ Σ, ker Df (p) intersects T p Σ transversally.Equivalently, d(det Df ), V = 0 on Σ. Definition 2. f is said to have a blowdown singularity along Σ if it only has corank one singularities and, in addition, ker Df (p) ⊂ T p Σ for every p ∈ Σ. Definition 3. f is said to have a cusp singularity along Σ if it only has corank one singularities and, in addition, at any point p ∈ Σ where d(det Df ), V = 0, i.e., where it fails to be a fold, one has d d (det Df, V ) , V = 0.

The case of a linear flight path
Without loss of generality, the trajectory of a straight flight path at height H > 0 and of constant, unit speed can be assumed to be γ(s) = (s, 0, H).For the straight flight path, R becomes R = (x 1 − s) 2 + x 2 2 + H 2 and the derivatives with respect to s become Proof of Theorem 2. We have Computing the determinant of the right projection π R : (x 1 , x 2 , s, τ ) → (x 1 , x 2 , ξ 1 , ξ 2 ), from ( 23) we obtain Thus, the differential Dπ R drops rank by one along the hypersurface s, and τ = 0 arbitrary }, and drops rank simply in the sense that d (det (Dπ R )) = 0 at Σ.We note that Σ is the set of points of C above the line Γ, the projection of the flight path on to the ground plane.In addition, with an isotropic antenna beam pattern, the entire problem is invariant with respect to reflection about the plane x 2 = 0, leading (as we will see) to a left-right artifact about Γ in reconstructions of the scene.
To classify the type of singularity of π R on Σ, one sees from ( 22), ( 23) that the kernel of Dπ R , necessarily spanned by a linear combination of the vector fields ∂ s , ∂ x 1 , ∂ x 2 and ∂ τ , has no ∂ x 1 , ∂ x 2 components and so is in fact spanned by a combination of only ∂ s and ∂ τ , both of which are tangent to Σ; thus, ker(Dπ R ) is tangent to Σ everywhere.Thus (see §5.2), π R has a blowdown singularity at Σ.
Similarly, for π L , from (14) one computes that the ∂ s and ∂ τ components of any vector in ker(Dπ L ) must be zero.The kernel is the null space in the 29) and (30) at x 2 = 0, we see that the second column of this matrix is zero, and thus ker(Dπ L ) = span{∂ x 2 }.This is transversal to Σ; hence, π L has a fold singularity at Σ (cf.§5.2).Thus, the canonical relation C is a fold/blowdown canonical relation as discussed in §4.This concludes the proof of Thm. 2.
To summarize: In the case of a straight flight path, the forward map F taking the scene V to the windowed DSAR data W is an FIO, given by ( 10), associated with a fold/blowdown canonical relation C.This means that, as discussed above, without beam forming to one side of the flight path or the other, backprojection will potentially create left-right artifacts in the image which are just as strong as the bona-fide part of the image.
Remark 2. As in Felea [11] it can be shown that an artifact relation, given by the graph of the canonical transformation χ(x, ξ) Here, I p,l (∆, Λ) is the class of pseudodifferential operators with singular symbols introduced in [26,22] now usually referred to as paired Lagrangian operators.As a result, an image extracted from F * F can have left-right artifacts just as strong as the actual features being imaged, i.e., the order of F * F is the same on ∆ and Λ.It is also possible to reduce the strength of these artifacts using a filtered backprojection method, with the principal symbol of the filter vanishing on Σ, along the lines of [17,33], but we will not pursue this here.
If, on the other hand, the system uses an antenna beam that illuminates only a region lying entirely to the left or to the right of the flight path, then F is an FIO associated with a canonical graph, which is a canonical relation satisfying the Bolker condition; consequently backprojection produces an image without artifacts.

The case of a circular flight path
The proof of Theorem 3 requires some preliminary discussion.7.1.Preliminaries.For simplicity, we consider the case of the flight trajectory being a circle of radius ρ > 0 at constant altitude in R 3 and centered above the origin in R 2 , parametrized by γ(s) = (ρ cos s, ρ sin s, ρh).Note that, in this case, we write the height H as H = ρh, where h is a dimensionless parameter (see Fig. 1).
We write e(s) = (cos s, sin s) and note that e ⊥ (s) := ė(s) = (− sin s, cos s).This yields From (11) we have, on the critical set, To determine whether artifacts can be avoided in the backprojected image, we consider the left projection, π L : C → T * R 2 .Let x and x correspond to the same s and ω.We consider the case when ω = ω and η Hence if there is an artifact for the circular trajectory case, that artifact must have R = R .(As it turns out, neither are they at points that are inverted with respect to the circle.)This is in contrast to standard monostatic SAR, since in that case R = R .Remark 3. The above shows that any artifacts arising from Doppler imaging for a circular flight path must appear in a location which is different from their location in standard monostatic SAR for the same flight path.Therefore, in principle, it might be possible to combine traditional monostatic SAR with Doppler SAR imaging in order to identify and remove artifacts.
In the next subsection, we investigate if it is possible to localize any backprojected artifacts to be outside the region of interest (ROI).7.2.Condition for C to be a local canonical graph.From the discussion in §3.1, we know that if C satisfies the Bolker condition (18) then backprojection results in an artifact-free image.We initially consider the first part of the Bolker condition, namely the requirement that C is a local canonical graph, and need to determine where the derivative of the left projection π L has full rank, i.e., rank(Dπ L ) = 4.We may parametrize C using coordinates (s, τ, x), with respect to which Next, to make the study of π L easier, introduce a new, s-dependent coordinate system in the plane, (u, v), defined by where Note that for s, x ranging over compact sets, u satisfies |u| ≤ 1 − for some > 0, and thus 1 − u 2 is bounded away from 0. The (u, v) coordinates are closely related to elliptic cylindrical coordinates, but in the plane and centered directly below the antenna; cf.[3, §2.7].We note that since e is perpendicular to e ⊥ , we have immediately (38) (x − ρe) • e = ρhS.
The new coordinate system is better understood once we have calculated Ṙ in the new coordinates, as follows.From ( 32) and (36) we have where where we have used C 2 − S 2 = 1.Hence, comparing (40) with (39), we find (41) Ṙ = −ρu.
Thus the coordinate u is simply proportional to the Doppler shift.Because γ = (ρe ⊥ , 0), from comparing ( 41) with ( 32), we see that u = R • (e ⊥ , 0), which is clearly bounded in magnitude by 1.The set u = ±1 corresponds to the line on the ground directly under the tangent to the flight path at e(s), and we exclude this by keeping the antenna beam pattern away from the direction of travel (forward and backward).
Observe that we may also parametrize C with the coordinates (s, τ, u, v).To see this, one needs to check that (42) clearly holds true since h > 0 means that both a and b in (43) are nonzero.
Therefore, by avoiding data from the forward and backward directions (or, alternatively, filtering out echoes associated to values of Ṙ near ±ρ), we have that (s, τ, u, v) forms a valid coordinate system on C. The coordinate system (s, τ, u, v) is designed to make any degeneracy of the projection π L appear in a single variable, namely v.In fact, in terms of (s, τ, u, v), one has To classify the singularities of π L , one needs to find first the set where Dπ L drops rank, i.e, where Computing R in the original coordinates x 1 , x 2 first, and then rewriting it in the new coordinates u, v using ( 29) and ( 30) one obtains Therefore, setting the v-derivative of (45) equal to zero shows that the set Σ of points where π L drops rank is given by where we have used (36).Since f (u 0 , v 0 ) = h 2 + 1 > 0, we see that the fiber of Σ lying over x = 0 is empty.
In particular, Dπ L is of full rank if and only if f = 0, where f is defined as in (46).Since Sh = (x−ρe)•e ρ , we have f > 0 iff and a sufficient, e-independent condition for this to hold is that Summarizing our analysis so far: Suppose that Σ is not in the microlocal support of F (for example, suppose that beam forming ensures that only points in D h,ρ are illuminated).Then the canonical relation of F is a local canonical graph.
Remark 5.If the flight path only consists of a portion of the full circle, then equation (51) below could be used to enlarge the region D h,ρ where the canonical relation of F is a local canonical graph.We do not pursue this idea further here.
Instead, in the next three subsections we analyze the structure of C near Σ.It will be convenient to introduce another set of coordinates as follows: Observe that (s, τ, p, q) also form a coordinate system on the canonical relation C. Indeed, since we have already remarked that (s, τ, x 1 , x 2 ) are coordinates on C, one only has to notice that for a fixed (s, τ )-value, the map (x 1 , x 2 ) → (p, q) is a diffeomorphism since the vectors d x p = e/ρ and d x q = e ⊥ /ρ are linearly independent.Points on Σ satisfy and hence is a defining equation for Σ.Since R 2 = ρ 2 (p 2 + q 2 + h 2 ), (51) can also be re-written as is also a defining equation for Σ.We will next make use of g to analyze the properties of Σ.
Proof of Theorem 3. Due to the length of the proof, it will be presented in parts: first we consider the singularity of π L (Sec.7.3), then the singularity of π R (Sec.7.4) and finally the elimination of artifacts (Sec.7.5).

7.3.
The singularity of π L .We showed in the last section that Dπ L drops rank at Σ.The following lemma establishes a key property of Σ.
Lemma 1.The function g in ( 52) is a defining function for Σ: dg = 0 at all points of Σ.
Proof.Clearly (∂g/∂q)(p, q) = 0 when q = 0, so we just need to check (∂g/∂p)(p, 0) = 0 whenever g(p, 0) = 0. To see this, we argue as follows: note that Recall that p = hS and suppose for contradiction that (∂g/∂p)(p, 0) = g(p, 0) = 0 for some p.Then Subsitituting ( 53) into (54) gives where we have cancelled a factor of hS, which is valid since (53) tells us that S = 0.The roots of this quadratic equation are Recalling that for points in Σ we have and substituting this into (55) gives which can only be true if we take the plus sign (since 0 < h, 0 < 1 − u 2 ).Simplifying (56), we get This leads to a contradiction since the left-hand side of (57) is strictly negative and it completes the proof of the Lemma.
Note that Lemma 1 shows that Σ is a smooth hypersurface in C.
If we illuminate a neighborhood of the boundary ∂D h,ρ , then it is of interest to know what kind of singularity π L has at Σ.This is answered by the next lemma.
Lemma 2. The projection π L has a fold singularity at Σ. Proof.Note that, at points of Σ, ker(Dπ L ) = span{∂ v } and we saw from the proof of Lemma 1 that dg = 0 there as well, so π L drops rank simply at Σ.Moreover, It follows that ker(Dπ L ) intersects T Σ transversally, proving the lemma.
Case 1: Assume β − = 0. Divide across (59) by hβ − to get where Let y := e v and multiply equation ( 60) by y to get the quadratic equation The solutions of this equation are From the definition of y, we must have y ± = e v ± > 0 for some value of v ± and therefore, y ± > 0.
Next we study the term under the square root in (61): Since η has the same sign as β − , the map V will be injective if since then, y + is the only possible positive root if β − < 0, while y − is the only possible root if β − > 0.
Lemma 5.The following inequality guarantees injectivity of V: Proof.In fact, if (64) holds, then from (38) we have Since 0 ≤ u 2 ≤ 1, this implies the inequality 1 − u 2 + 2hS < h 2 ; the left of which is decreased by multiplying by 1−u 2 to obtain Moving h 2 S 2 to the left side, we obtain (1−u 2 +hS) 2 < h 2 C 2 , which can be rewritten . This is exactly condition (62), as claimed, so that (64) is sufficient for injectivity of V.
Writing (64) as x • e < ρ + ρ(h 2 − 1)/2 and considering all possible locations on the flight path, it then follows that a sufficient condition to guarantee injectivity of V is This concludes the proof of Theorem 3. Remark 6.Note that condition (66) implies that the left projection π L is injective, and also guarantees that condition (47) is automatically satisfied, so that π L is an immersion.Therefore, C is a canonical graph over the region defined by (66), ensuring artifact-free imaging of scenes there via filtered backprojection.
Remark 7. To summarize the microlocal analysis so far of DSAR for a circular flight path, C is a canonical relation whose projections, π L and π R , are of fold and cusp type, resp.A similar geometry appeared in [16] in the case of monostatic SAR when the flight path had simple inflection points.In that situation it was shown that F * F produces an artifact relation which is a canonical relation, with a codimension-two set of points where it is nonsmooth, called an open umbrella; see [15] for background material on umbrellas and their relevance in seismic imaging, and [16] for how this geometry arises in monostatic SAR.

First order correction to the start-stop approximation
We now describe and analyze a correction to the start-stop approximation, Assumption 3, which in the form of R(t + T tot ) = R(t) was substituted into (4) and used to derive (7) and thus (10).The correction we now consider comes from modifying the discussion below (3) by including a first order term in the expansion of R(t + T tot ).See also [39,7] for discussions of the start-stop approximation and its limitations.
As in Sec. 2, T tot is determined implicitly by ( 4), namely c 0 T tot = R(t)+R (t + T tot ).However, we now expand Solving for T tot and ignoring terms O(c −3 0 ), we obtain the first order refinement of the start-stop approximation, namely that the total travel time is (67) and thus, under this refined approximation, the scattered wave arrives at time Substituting into (68) the linear terms of the Taylor expansions and, calculating modulo (t − s) 2 , yields the approximation W 1 of W analogous to (10): (70) W 1 (s, ω) = e iφ(s,ω,x;τ ) a(x, s; τ )V (x)dτ dx, where Rewriting the expression in (71) that multiplies c −2 0 as (1/2) (R 2 ) .. for compactness (where the superscript upper right dots still denote differentiation in t), one sees that this modified phase function φ parametrizes the canonical relation The fold/blowdown structure of the canonical relation in the case of the linear flight track, established in §6, is not structurally stable, since arbitrarily small perturbations of blowdowns are not (in general) blowdowns.Nevertheless, we now show that under the first order correction to the start-stop approximation, C mod is still a fold/blowdown, indicating the robustness of our approach.
To see this, note that the analogue of ( 23) taking the correction into account is which vanishes to first order at (s, x, τ ) ∈ Σ = {x 2 = 0}.Furthermore, at these points, ker(Dπ R ) = R • ∂ ∂s ⊂ T Σ, so that π R has a blowdown singularity at Σ.A similar calculation shows that Dπ L has a kernel whose nonzero elements have a nonzero coefficient of ∂ ∂x 2 , so that π L has a fold singularity at Σ, and C mod is a fold/blowdown canonical relation.
In the case of the circular flight track, the fold/cusp structure of the canonical relation derived in §7 under the start-stop approximation, is structurally stable in the following sense: a small perturbation of the phase function (say in the C 4 topology) will result in a small perturbation of C and this causes small C 3 perturbations of the projections π L and π R .Folds are stable under C 2 perturbations and cusps are stable under C 3 perturbations [18], and thus one expects that the first order correction to the start-stop approximation will not essentially change the microlocal analysis for a circular flight path: C mod will still be a fold/cusp, so that the artifacts are of the same type and strength as shown above, although with their locations moved slightly.However, we have not proven this and it is a subject for future research.

Concluding Remarks
We briefly compare our results for Doppler SAR with the case of monostatic SAR treated in [32].In all cases, we can expect the strength of any associated artifacts in the image to be as strong as the bona-fide part.
For the case of a linear flight path, the results for Doppler are the same as for monostatic SAR, in the sense that π L and π R have fold and blowdown singularities, resp.The normal operator F * F formed without beam forming will have the same strong artifact, i.e., F * F ∈ I −1,0 (∆, C χ ), with C χ a canonical graph (cf.Remark 2 in Sec. 6).
On the other hand, for a circular flight path, the microlocal geometry for Doppler SAR differs from that for monostatic SAR: we have shown that in the case of Doppler, the singularities of π L , π R , are of fold and cusp type, resp., whereas for monstatic SAR, both singularities are folds [32].Regarding the operator F * F, it was shown in [11] that in the monostatic SAR, F * F ∈ I −1,0 (∆, C), where C is a two-sided fold.For Doppler SAR we expect, based on the results of Felea and Nolan [16], that F * F ∈ I −1 (∆, C), where C is an open umbrella (see Remark 7).For the circular flight path, we have described an explicit region D h,ρ where the wave front relation of F is a canonical graph, so that no artifacts appear at the back projection, allowing accurate imaging of the terrain.In addition, as described at the end of Sec.7.1, the artifacts arising from Doppler imaging in a circular flight path geometry are spatially separated from their location for monostatic SAR for the same flight path.It should therefore be possible to identify and remove artifacts by combining monostatic SAR data with Doppler SAR data; we hope to return to this in the future.

Figure 1 .
Figure 1.Schematic of a circular flight path.

Figure 2 .
Figure 2. The curves of constant u (hyperbolas) and constant v (vertical lines) for a location on the flight path in which the flight velocity vector is along the vertical axis.The coordinate system is centered directly under the antenna.