The Hyperbolic Geometry of SAR Imaging

Andrew S. Milman

Unfortunately, the Netscape does not display the Greek letters correctly; to see them, you will have to use Internet Explorer.

July 1, 1999

Abstract

This paper shows how the use of hyperbolic functions allows us to write an exact mathematical representation of the SAR imaging. This problem is primarily a geometric one, that of accounting for curved wavefronts: the spirit of this paper is to emphasize these geometrical properties over electromagnetic ones.  This gives us a new and fruitful way to think about SAR imaging.  Within this framework, I show how to correct for deviations of a radar from a straight flight path. This method will work even in situations where the curvature of the wavefronts is very large, where traditional methods do not. The image-formation algorithm, called w-k migration, that results from this analysis of SAR imaging is simpler and faster than polar formatting, especially for radars with very large beamwidths. As an added benefit, w-k migration is surprisingly simple to derive.

1. Introduction

Historically, processing synthetic aperture radar (SAR) data to produce images has used an approximation to express the instantaneous distance from the radar to a point on the ground. While this approximation is often adequate, it is still an approximation: we would prefer to have an exact expression both to aid our understanding of the SAR imaging problem and to improve the accuracy of our computations. In this paper, I derive an exact mathematical representation of the SAR-imaging problem. We should be pleased to see that the result is very simple. In particular, I shall show how, within this framework, we can make exact corrections for lateral motions of the radar.

      The usual SAR-imaging approximation ignores the curvature of the wavefronts: this leads directly to the now-familiar polar formatting (see Walker, 1980; Jakowatz et al. provide a good discussion of SAR processing using polar formatting). When the SAR antenna has a narrow beamwidth, this approximation may be adequate, but it breaks down when the beamwidth is large, as it will be for low-frequency SAR’s of the type that are being studied for foliage-penetrating applications. However, using hyperbolic geometry, we can develop expressions that take the curvature into account exactly:  no approximation is necessary.  This leads to an alternative method of producing SAR images, which is called w-k migration (see Milman, 1993, and references cited therein; Prati et al., 1990; and Cafforio et al., 1989, 1991a and 1991b, are the original papers on the subject, as far as I know; see also Soumekh, 1990). Not only is this method more accurate than polar formatting, it is also simpler and easier to implement on a computer.  There are many advantages and no known drawbacks to using w-k migration.

      When we have a SAR with a very wide beamwidth, traditional methods of compensating for lateral motions of the radar.  Normally, we compensate for these motions of the SAR by multiplying each pulse by the appropriate complex exponential. This works when the curvature is small—the antenna beam is narrow—but it doesn’t work for wide beams. I use the w-k migration formalism to develop an exact correction that works no matter how large the curvature is.

      Polar formatting has the limitation that, because of the plane-wave approximation, the focus of an image deteriorates with increasing distance from the scene center. This requires that we form many sub-images and patch them together (see Figure 1). Aside from the computation costs, this introduces discontinuities in the phase of the complex image; these discontinuities will affect algorithms that use this complex phase. We shall see that w-k migration processing has none of these drawbacks.

      One way to view migration processing is to consider, as I do here, the geometrical aspects of the problem, rather than concentrating on the electromagnetic formalism.  If we treat electromagnetic waves as being planar, we can associate electromagnetic frequency with range, at least for a chirped radar. This concept leads to the usual polar-formatting diagram, where pulses are laid out in a radial manner, with the near and far edges of the swath being taken as corresponding to the frequencies f0 and f0 + B, the minimum and maximum frequencies in the chirp (i.e., B is the bandwidth of the radar; see Figure 2).  The curvature of the lines of constant range, in this picture, depends on the ratio of the bandwidth B to the center frequency f0 + ½B. 

      In migration processing, we focus on the geometrical aspects of the problem: we shall consider a pulse as providing information about the backscattered intensity as a function of range, rather than frequency, with the range having an absolute meaning—distance from the radar to a target; the frequency of the pulse plays only a minor role.  In the development that follows, I take the fundamental measurements as being the amplitude and phase of the backscattered signal as a function of pulse number—or aircraft position—and range. The first step in the analysis requires us to perform a 2D FFT of the data; let k be the range wavenumber. If the ranges to the near and far edges of the swath are r0 and r1, with a resolution dr = c/(2B)—where c is the speed of light—then the range wavenumber k lies in the interval (0, 1/dr), not (f0, f0 + B). Transmitted wavefronts are spheres centered on the radar; we can account for the curvature exactly if we use the right coördinate system. This may be a subtle point, but it is central to understanding the difference between the method I am proposing here and work by other authors (this geometric viewpoint is inherent in Milman, 1993).

      Carrara et al. (1995, pp. 401 ff) and Golden et al. (1994) discuss a similar method that they call range migration.  Although it is based, in part, on Milman (1993), they do not make use of the property that the migration formalism is an exact representation of the SAR imaging problem. This may be, in part, because they have not adopted the geometrical outlook I espouse here.

      There is always a need to correct for lateral motions of the radar, i.e., known deviations from a straight flight path. I shall call such corrections motion compensation in this paper, although this term is not universally recognized in the SAR community.  In this paper, I show how to correct for these motions globally, in the sense that, assuming that we know the deviations from the nominal, straight flight path, we can perform a single correction to the data that compensates for that motion everywhere within the field-of-view of the radar.  By contrast, the “range migration” algorithm discussed by Carrara et al. (1995) and Golden et al. (1994) (also personal communication) compensates for lateral motions by creating sub-images with the data compensated to the center of the scene.  This only accounts for the lateral motions of the radar approximately and requires that the same raw data be processed several times, once for each sub-image that the signal history contributes to.

      One attractive aspect of using w-k migration is it gives us a new perspective on the problem and a new way of thinking about it.  In this paper, I show how we can make use of hyperbolic functions in the context of correcting for deviations from a straight flight path.

Figure 1.  A) Forming several subimages; B) forming a single image.

2. An overview of w-k migration

Before I start the derivation of w-k migration, I should show you where we’re going.  I want to show that there is a closed-form expression that relates the signal history from the SAR to the distribution of targets on the ground, and that this expression has an analytical inverse.  In  so doing, I make use of Hankel transforms, which are similar in many ways to the more familiar Fourier transform.  I discuss Hankel transforms in more detail in Appendix A. 

      We can write a Fourier transform as

                                                  ,                                            

where g(x) is an arbitrary function and the Fourier kernel KF = exp{‑2piwx}.  This transform has an inverse

                                                 ,                                            

where KF* is the complex conjugate of KF.  Similarly, a Hankel transform is                              

                                                  ,                                            

where the Hankel kernel is

                                                         ,                                                    

and where H0(2)(x) = J0(x) – iY0(x) is a Hankel function of zero order of the second kind; J0 and Y0 are Bessel functions of the first and second kind. This transform also has an inverse:

                                                 .                                            

Hankel transforms are less familiar than Fourier transforms but have similar properties.

      Now suppose I have a function A(x) whose Fourier transform is a(w).  If I measure a(w), it is trivial to calculate A(x): we just use an inverse FFT.  Similarly, if B(x) is the Hankel transform of b(w), we can calculate B from b with an inverse Hankel transform.  We shall see later that we can use the usual FFT software to calculate a Hankel transform, so it is as efficient to calculate as a Fourier transform.

      Suppose that s(t,r) is the measured signal history from a SAR, where r is distance from the radar to a reflector and t, the (slow) time that measures the position of the SAR along the flight path (I am using units such that the speed of the SAR, v = 1).  I shall define S(w,k) to be the 2D Fourier transform of s(t,r) and show that S(w,ky) is a 2D transform of the scene reflectivity s(x,y), where ky2 = k2 - w2. The change of the independent variable from k to ky is an important part of this theory.  The desired transform is a Fourier transform in x, the along-track direction, and a Hankel transform in y, the cross-track direction.  Let me write the combined Fourier-Hankel kernel as

                                               .                                          

Then I shall show that

                                   .                              

But if this is so, we can form a SAR image easily: the image is just

                                  .                             

      This is an exact relationship: I have not had to make any plane-wave or similar approximation to arrive at equation .  Furthermore, this relationship is valid regardless of the width of the size of the region illuminated by the antenna; there is no loss of focus at the edges.  I am always surprised at how simple this solution is.

      It is customary to analyze many kinds of SAR imaging problems by writing the signal history as a Taylor series and keep the first two or three terms.  This is often an analytical nightmare, and it requires us to justify ignoring all terms higher than a certain order. But because w-k migration provides an exact analytical relationship between signal history S(w,ky) and the scene s(x,y), we can analyze SAR imaging problems quite easily within this framework.  I provide one such analysis in this paper, pertaining to the question of how to correct for lateral motions of the SAR, but first, I shall provide a rigorous derivation of the results that I have claimed in this section.

3. Geometry and signal history for a straight flight path

Consider a radar mounted on an airplane that is used to make synthetic aperture radar (SAR) images of the earth’s surface. Something that is well known, but apparently often overlooked, is that the shape of the curve that describes the instantaneous distance from the SAR to a fixed point on the ground is an hyperbola. We can make use of this property to simplify the analysis of the SAR imaging problem; Figure 3 shows the relevant geometry.  Here, h is the height of the radar above the (plane) surface of the earth.  The radar moves with velocity v along the x-axis. Assume that the time t=0 corresponds to x=0, so vt measures the position of the SAR along the flight track.  To make the notation simpler, I shall use units chosen so that v = 1; therefore, t is a distance (as well as slow time). The distance from the SAR, located at the position [t, 0, h], to a point at (x, y, 0) on the surface is

                                                    .                                              

We can reduce this to an equivalent 2-dimensional problem with the distance in the ground plane given by

                                               .                                       

Then the signal history is

                                                 ,                                         

where s0(x,y) is the complex radar backscattering coefficient at location (x,y) and the curve C is the path defined by r2 = y2 + (x ‑ t)2.

      The radar wavelength is l; I assume that it has been modified appropriately to account for the projection into the ground plane.  Alternatively, we can take this to be the signal history that would have been received by a hypothetical SAR flying at zero altitude. It makes no difference for the purposes of this discussion. The complex exponential represents the relative phase of the reflected signal.

      The path of a point on the ground, in the frame of reference fixed to the airplane, is

                                                       ,                                               

which, of course, is the equation for an hyperbola.  This suggests that we use hyperbolic functions to describe the geometry of SAR imaging. I shall show how this allows us to use an exact analytic description of the SAR geometry.  This will lead to an algorithm for producing SAR images that is mathematically exact; this is, in part, a simplified derivation of the so-called w-k migration algorithm described by Milman (1993) and others.

Figure 2.  Polar formatting, showing the original data (solid lines), which form an angle f corresponding to the angle subtended by the synthesized aperture, and f0 and B, which are the lowest frequency in the chirped pulse and the bandwidth. Each pulse is plotted in a radial direction. The dotted lines are the data interpolated to the x-y coördinate system.

Figure 3. SAR imaging geometry.

      The hyperbolic functions sinhy and coshy are related by

                                                     .                                             

If we let

                                                           ,                                                   

then

                             .                     

Figure 4 shows the geometric meaning of the angle y. It shows the path of an object on the ground located at a distance y from the SAR track, as seen from the aircraft. Now, since y =rsechy,  we can write x and y in the form

                                                          ,                                                  

                                                             .                                                     

      We can develop an exact representation the phase history using hyperbolic functions.  We will need to transform the variables of integration from (x,y) to (r,y).  Using the relations tanh2y + sech2y = 1,

                     ,              

we see that the Jacobian of this transformation is rsechy.  I shall write the signal history as

                    .            

The factor 1/r was introduced to cancel the factor r in the Jacobian.

      Now take the 2-dimensional Fourier transform of s(t,r):

                                          

                    .             

Transform this by letting y = rsechy, so

                  .          

Now substitute x = t + ysinhy, so

Figure 4.  Hyperbolic coördinates.

                  

                             .                     

4. Stolt and Hankel transforms

We now can make the following substitutions: let

                                                      ,                                               

                                                            ,                                                    

and

                                                             .                                                   

Remember that k is the wavenumber in the radial direction; this substitution will make the transformation to ky, the wavenumber in the y-direction. By doing this, we effectively remove the curvature of the wavefronts. Using the relationship coshycoshf + sinhysinhf = cosh(y+f), we get the simple result that

                            .                     

The transformation of S(w,k) to S(w,ky) is a one-dimensional interpolation, called the Stolt transform.

      Compare equation , which is the integral representation of a Hankel function to the integral over t in equation ;  this shows that equation is a Hankel transform of the complex reflectivity s0. I describe the relevant properties of Hankel functions and Hankel transforms in Appendix A.

      Knowing this, we can perform the integration over t analytically, using equation :

                                 

                                       .                               

The inverse of equation is

                           .                   

      We should note the factor exp{-2piykycosh(y‑f)} in equation :  to the extent that we can take the infinite limits of integration seriously, the value of f is immaterial.  In reality, there is only a finite angle y that we integrate over, so the value of f will make a small difference. However, this is a limitation on our ability to form images from a finite amount of data, and will arise no matter how we process the data: it is not a shortcoming of the w-k migration method.

      This shows that, without making any plane-wave approximation, it is possible to write the radar cross-section distribution in terms of the signal history. We can invert it by using a Fourier transform in the along-track direction, and a Hankel transform in the cross-track direction.  The Hankel transform is unfamiliar and (to my knowledge), no one has Fast Hankel Transform (FHT) software the way we have FFT’s.  Fortunately, we can use the approximation given by equation in equation , which is accurate as long as 2pyky > 3 (see Appendix A). So the Hankel transform we need to perform is, aside from division by (2pyky)½, equivalent to an FFT, and we can use the usual FFT software.

      The form of the Stolt transform given in equation is different from the usual definition—e.g., Milman (1993).  Here, I have used k + 2/l, instead of just k, to account formally for the factor exp{‑4pir/l} in equations et seq. There are different ways that we might take this into account in practical problems, depending on the specifics of the situation: the difference between using k and k + 2/l is formal, rather than substantive.

5. Correction for lateral motions

Figure 5.  If the SAR is displaced by an amount q from the x-axis to the x¢-axis, we should correct the measurements from the line c¢‑d¢ to c‑d.  Changing the phase of each pulse shifts the data radially from r¢ to r.

Now what happens if the SAR drifts laterally from a straight flight line? Suppose that, at time t, the SAR is displaced by -q(t) in the y-direction.  This is equivalent of having the scene move in the y-direction by +q(t)—see Figure 5.  If we simply multiplied each pulse by a complex phase to correct for these motions by shifting it toward or away the SAR by an amount ‑q, we would be making an error because, while the SAR was displaced in the y-direction, the correction is being applied in the radial direction. For wide-beamwidth SAR’s, this is a serious problem. We can perform the correction properly as follows.

      When we include the lateral motions, equation becomes

                         .                  

We correct for the later motion analytically as follows. Remember that we will define ky = [(k + 2/l)2‑w2]½; define S(w,k) to be the Fourier transform of exp{2pikyq(t)} s(t,r): 

                                       

                                                          

[compare this with equation ]. This complex factor will correct for the lateral motions denoted by q(t). Note that, because of the factor ky in this expression, rather than k, we are making a correction in the y-direction, and not the radial direction as the normal motion-compensation would have done.  This step is crucial when the beamwidth is large: we would be making a large error if we compensated for lateral motions of the SAR—which by definition are in the y-direction—by making corrections in the radial direction. To continue, substitute y = rsechy, and then x = t + rsinhy:

           

[compare this with equation ] where I have used the notation q¢(x) º  q(x +rsinhy) for simplicity. Now substitute y¢ = y ‑ q¢(x):

          .   

Making the same substitutions as before for k and w [equation ],                                        

                .        

      Now use equation to perform the integration over y:

                  ,           

and finally, use equation to approximate H0(2):

                    .            

Since q is a few meters and r is tens of kilometers, we can safely replace [r‑q¢(t¢)]‑½ with r-½ in this equation; it affects the amplitude of the signal, but not the phase. Finally, we see that the two exponentials that contain q¢(x) cancel—which is why I introduced the factor exp{2pikyq(x)} into equation in the first place.  So

                               .                        

Therefore, we have corrected for the lateral motions of the SAR in a way that works just as well when the wavefront curvature is large as when it is small. Note that I have made no approximations to the geometry of the situation, I only made an approximation to the Hankel function to speed up the data processing.

6. Autofocus

A problem similar to the one I discussed in the last section can arise with autofocus, which involves correcting for the phase change induced by lateral motions that are smaller than a resolution cell (e.g., see Jakowatz et al., 1996).  Unlike the correction for the larger-scale lateral motions, which apply a constant phase rate to each pulse to use the Fourier shift theorem to correct for the motion, this is a single constant phase shift that applies to all the data within a single pulse.  Since it varies from pulse to pulse, it affects the azimuth focus of the data. Most autofocus algorithms estimate the appropriate phase shift for each pulse from the image data and apply the correction in the y-direction. However, because of the geometrical transformation inherent in polar formatting (see Figure 2), the original range direction is no longer parallel to the y-direction in the polar-formatted data. Therefore the autofocus correction, if applied to the polar-formatted data, will apply the corrections to the wrong pulses: this error is small near the center of the image, but grows with increasing distance from the center.  This limits the accuracy of the autofocus algorithm and could even prevent it from converging.  Since w-k migration processing maps each pulse to a line in the y-direction, this error does not arise.

7. Conclusion

There is a straightforward method for performing the motion compensation for SAR imaging. It can also be used in conjunction with autofocus techniques to focus the images.   Like other aspects of w-k migration processing, it is based on an exact mathematical relationship.  It also shows how viewing the SAR imaging process in terms of hyperbolic functions simplifies both the basic image-formation algorithm and the corrections to the data necessitated by deviations from a perfectly straight flight path. In general, the w‑k migration algorithm is simpler than polar formatting and does not have the limitations on image size that are imposed in polar formatting by the plane-wave approximation.

      Although I have derived the motion-compensation algorithm here only for horizontal deviations from a straight flight path, this method can also be applied to the problem of out-of-plane motions.

Appendix A.  Hankel Functions

We shall use the zero-order Hankel functions, which are defined as H0(1)(z) = J0(z) + iY0(z) and H0(2)(z) = J0(z) ‑ iY0(z), where J0 and Y0 are the zero-order Bessel functions of the first and second kind. We shall use the integral representations

                                                    ,                                             

                                                   .                                            

      Later, we will use the following asymptotic approximations to H0(1) and H0(2):

                                                    ,

                                                   .                                            

      These approximations for H0(1)(x) and H0(2)(x) are very accurate for x>3. We shall see that this will always obtain for any plausible SAR. In the event that we need a more accurate asymptotic expansion, a more elaborate one is given below.

      We can write the inverse Hankel transform as follows.  If F(r) is a function such that

                                                         ,                                                 

then

                                     .                              

(This relationship is also true when H0(1)(ur)H0(2)(us) is replaced with C0(ur)C0(us) and C0 is either of the Bessel functions J0 or Y0.). This is just a restatement of equations and .

      A more accurate asymptotic expansion is the following: 

                       ,               

where n is the order of the Hankel function and

                                  

and