Synthetic aperture radar image compression

Abstract

A method is disclosed of generating multi-resolution images using synthetic aperture radar (SAR) phase history data. A high resolution SAR image may be formed from the phase history data, or a lower resolution SAR image may be formed from a sub-aperture of the phase history data. The multi-resolution images may also be compressed. Compressed multi-resolution images are progressively transmitted to a client until the client receives an image with adequate image quality and/or resolution. Compressed multi-resolution images may also be used for iterative target detection, where images are analyzed by a target detection processor, starting with a lowest resolution image, going to a highest resolution image, until an adequate image of the target or scene is obtained. If a target is not found in the multi-resolution images, a new sub-aperture is chosen on the phase history data, and the target detection process repeats for the new sub-aperture.

Description

TECHNICAL FIELD

The present invention relates generally to processing data from a synthetic aperture radar (SAR) imaging system and, more particularly, to a system and method for generating compressed, multi-resolution SAR images.

BACKGROUND OF THE INVENTION

Synthetic aperture radar (SAR) is a well known technique for developing radar imagery with excellent two-dimensional (2D) resolution. This is typically done by flying an airplane over the ground to be mapped, and successively transmitting a sequence of radar pulses. With the forward motion of the airplane, each successive radar pulse is transmitted from a position a little further along on the flight path, thus simulating a very long radar array.

Synthetic aperture radar is now a mature technique used to generate radar images in which fine detail can be resolved, and provides unique capabilities as an imaging tool. Because it provides its own illumination (the radar pulses), it can image at any time of day or night, regardless of sunlight illumination. And because the radar wavelengths are much longer than those of visible or infrared light, SAR can also “see” through cloudy and dusty conditions that visible and infrared instruments cannot.

As the radar moves, a pulse is transmitted at each position; the return echoes pass through the receiver and are recorded in an “echo store.” Because the radar is moving relative to the ground, the returned echoes are Doppler-shifted (negatively as the radar approaches a target; positively as it moves away). Comparing the Doppler-shifted frequencies to a reference frequency allows many returned signals to be “focused” on a single point, effectively increasing the length of the antenna that is imaging that particular point. This focusing operation, commonly known as SAR processing, is now done digitally on powerful computer systems. SAR processing must correctly match the variation in Doppler frequency for each point in the image, however, and this requires very precise knowledge of the relative motion between the platform and the imaged objects (which is the cause of the Doppler variation in the first place).

Radar energy is transmitted in the form of sequential pulses, at different time instances in the flight path of the vehicle. The pulses interact with the terrain (and any object on the terrain) and a portion of the pulse energy is reflected back towards the platform and recorded by a detector. Returns from different objects arrive at different times at the detector. These time differences provide information on range, which is then used to create a final radar image.

Typically, the transmitted pulses in a SAR system are sinusoidal in nature. The detector records the reflected sinusoids by recording a complex number, whose phase component is directly related to the time-of-flight, or range, and its magnitude component is proportional to the reflected energy. Thus, data captured in a SAR system is complex in nature, with properties that are unique with respect to data captured in systems from other imaging modalities.

FIG. 1 is an exemplary SAR image capture system according to the prior art. In block 10, a pulse generator generates chirps to be emitted from the SAR antenna. Before the chirps are emitted by a transmitter in block 12, however, they are frequency modulated in block 11. A receiver receives chirp echoes (i.e., reflections of the frequency modulated chirps from objects in the area being imaged), and block 13 subsequently performs SAR phase history data preprocessing. The data is then digitized in block 14. Block 15 may then reconstruct the two-dimensional SAR phase history data by performing such steps as pulse compression, Doppler shift corrections, and the like. A two-dimensional (2D) SAR image may then be obtained in block 16 by taking the inverse Fourier transform of the reconstructed phase history data.

Generally, return signals from the transmitted pulses are sampled in the airplane and either processed on board for immediate exploitation or stored or transmitted for processing at another site. The processing is computationally expensive, employing such techniques as FFT, inverse FFT, or correlation on vast amounts of data. These operations require vast processing power and storage. FIG. 2A, for example, illustrates the formation of SAR image 206 by taking the inverse FFT (IFFT) 204 of SAR phase history data 202. Consequently, SAR phase history data 202 may be obtained by taking the FFT 208 of SAR image 206. Furthermore, a compressed SAR image 216 may be obtained by taking the IFFT 214 of compressed SAR phase history data 212, whereby compressed SAR phase history data 212 may be obtained by taking the FFT 218 of compressed SAR image 216. Current methods of SAR phase history data compression are discussed in A Review of Current Raw SAR Data Compression Techniques by A. El Boustuni et al.

Technology trends in the field of SAR indicate that SAR system designs are continuously pushing the envelope for increases in area coverage and resolution. These trends imply massive amounts of collected data, which in turn, stress the ability to store collected SAR data and rapidly disseminate SAR data. Furthermore, SAR data transmission is undesirably time-consuming for customers who are only interested in a particular target in the scene, as opposed to the entire scene, or customers who aren't interested in a high resolution image of the scene.

SUMMARY OF THE INVENTION

The present invention is embodied in a method of generating multi-resolution synthetic aperture radar (SAR) images. The method includes the steps of: (a) obtaining SAR phase history data, denoted as PH; (b) calculating an inverse Fourier transform of a first data set, denoted as ph₀, of the SAR phase history data, PH, to form a highest resolution image, denoted as G₀, where G₀=FT⁻¹(ph₀); (c) calculating the inverse Fourier transform of a plurality of other data sets, denoted as ph₁ . . . ph_N, to form a plurality of correlated lower resolution images, denoted as G₁ . . . G_N, where G_i=FT⁻¹(ph_i), for i=1 to N; (d) subtracting successive ones of the highest resolution image, G₀, and the lower resolution images, G₁ . . . G_N, to obtain a plurality of edge enhanced images, denoted as L₀ . . . L_N−1, where L_i=G_i−G_i+1, for i=0 to N−1; and (e) obtaining a lowest resolution image, denoted as L_N.

In a further embodiment, the present invention includes a method of generating multi-resolution synthetic aperture radar (SAR) images. The method includes the steps of: (a) obtaining SAR phase history data, denoted as PH; (b) computing an inverse Fourier transform of a first data set, denoted as ph₀, of the SAR phase history data, PH, to form a highest resolution image, denoted as G₀, where G₀=FT⁻¹(ph₀); (c) successively low-pass filtering the highest resolution image, G₀, to form a plurality of correlated lower resolution images, denoted as G₁ . . . G_N, where G_i+1=LPF(G_i), for i=0 to N−1; (d) subtracting successive ones of the highest resolution image, G₀, and the lower resolution images, G₁ . . . G_N, to obtain a plurality of edge enhanced images, denoted as L₀ . . . L_N−1, where L_i=G_i−G_i+1, for i=0 to N−1; and (e) obtaining a lowest resolution image, denoted as L_N.

In yet another embodiment, the present invention includes a method of iterative target detection in a synthetic aperture radar (SAR) scene. The method includes (a) obtaining SAR phase history data, denoted as PH, corresponding to the scene; (b) integrating over a plurality of data sets, denoted as ph₀ . . . ph_N, of the phase history data, PH, to obtain a plurality of features, denoted as F₀ . . . F_N; (c) storing the plurality of features, F₀ . . . F_N; (d) performing target detection based on the stored plurality of features to obtain a receiver operating characteristic (ROC) of the scene; (e) determining if the ROC is acceptable; and (f) repeating steps (b)–(e) for a new plurality of other data sets if the ROC is unacceptable. In still another embodiment, the present invention includes a method of progressively transmitting multi-resolution synthetic aperture radar (SAR) images. The method includes the steps of: (a) obtaining SAR phase history data corresponding to a scene; (b) generating a correlated image pyramid including a plurality of multi-resolution images of the scene; (c) generating a Laplacian image pyramid including a plurality of Laplacian images corresponding to a difference of the multi-resolution images in successive levels of the correlated image pyramid; (d) encoding each one of the plurality of Laplacian images; (e) receiving, from a client, a request for a desired image of the scene; and (f) transmitting each one of the plurality of Laplacian images, starting from a bottom level of the Laplacian image pyramid, until the desired image of the scene is transmitted.

It is to be understood that both the foregoing general description and the following detailed description are exemplary, but are not restrictive, of the invention.

BRIEF DESCRIPTION OF THE DRAWING

The invention is best understood from the following detailed description when read in connection with the accompanying drawing. It is emphasized that, according to common practice, the various features of the drawing are not to scale. On the contrary, the dimensions of the various features are arbitrarily expanded or reduced for clarity. Included in the drawing are the following figures:

FIG. 1 is a flow diagram of a SAR image acquisition system, according to the prior art;

FIG. 2A is a block diagram of an exemplary SAR image generator, according to the prior art;

FIG. 2B is a block diagram of an exemplary compressed SAR image generator;

FIG. 3 is a pictorial depicting a comparison between forming a high resolution image using an entire phase history data set and forming a low resolution image from a sub-aperture on the phase history data set;

FIG. 4A is a plot of SAR phase history data with exemplary rectangular sub-apertures, according to the present invention;

FIG. 4B is a plot of an exemplary three-dimensional RECT sub-aperture function, according to the present invention;

FIG. 5A is a plot of SAR phase history data with exemplary circular sub-apertures, according to the present invention;

FIG. 5B is a plot of an exemplary three-dimensional Gaussian sub-aperture function, according to the present invention;

FIG. 6 is a plot of SAR phase history data with exemplary sectoral sub-apertures, according to the present invention;

FIG. 7A is a pictorial depicting the formation of a one-dimensional Gaussian pyramid for complex data points, according to the present invention;

FIG. 7B is a plot of exemplary Gaussian weighting functions for use in obtaining the “g₁” level in the one-dimensional Gaussian pyramid of FIG. 11A, according to the present invention;

FIG. 7C is a plot of another exemplary Gaussian weighting function for use in obtaining the “g₂” level in the one-dimensional Gaussian pyramid of FIG. 11A, according to the present invention;

FIG. 8 is a block diagram of an exemplary embodiment of a system for Laplacian pyramid coding and decoding of multi-resolution SAR images, according to the present invention;

FIG. 9 is a flow diagram illustrating a method of generating compressed multi-resolution SAR images, according to the present invention;

FIG. 10 is a pictorial depicting an exemplary system for generating compressed multi-resolution SAR images from a high resolution image using all the phase history data and a low resolution image from sub-aperture phase history data, according to the present invention;

FIG. 11 is a pictorial depicting a system for generating compressed multi-resolution SAR images from circular sub-apertures, according to the present invention;

FIG. 12 is a pictorial depicting a system for generating compressed multi-resolution SAR images from rectangular sub-apertures, according to the present invention;

FIG. 13 is a pictorial illustrating an exemplary target detection system using sectoral sub-apertures, according to the present invention.

DETAILED DESCRIPTION OF THE INVENTION

Airborne and spaceborne based SAR imaging systems are active coherent image acquisition systems and work like flashlight photography. A SAR imaging system emits a microwave frequency-modulated chirp and receives its echo. The collected raw data are then digitally processed to generate 2D SAR phase history data and a SAR image.

More specifically, a SAR system consists of the following subsystems: a pulse generator for generating particular microwave wavelength pulses, an antenna for transmitting the generated pulse to illuminate a region of interest, an antenna for receiving the backscattered returned pulses (i.e., echoes) from target points in the region of interest, a digitizer for converting the received complex analog signal (amplitude and phase) to a digital format, and a signal processor for constructing the SAR 2D phase history data and the SAR image.

SAR signal preprocessing and postprocessing are based on convolution of the SAR emitted beam with the scene and de-convolution of the collected data into a high quality 2D image. To describe the process of constructing the phase history data and forming the SAR image, a few relevant terms are defined, below.

Discrete Fourier Transform (DFT):

$\mathrm{DFT}\left[h\left(m,n\right)\right]=H\left({\omega }_k,{\omega }_l\right)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}{e}^{-i\left({\omega }_{k}m+{\omega }_{l}n\right)}h\left(m,n\right)$ ${\mathrm{DFT}}^{-1}\left[H\left({\omega }_k,{\omega }_l\right)\right]=h\left(m,n\right)=\frac{1}{\mathrm{MN}}\sum _{k=0}^{M-1}\sum _{l=0}^{N-1}{e}^{+i\left({\omega }_{k}m+{\omega }_{l}n\right)}H\left({\omega }_k,{\omega }_l\right)$

Those skilled in the art will recognize that the DFT may be efficiently implemented with the Fast Fourier Transform (FFT).

The convolution of two signals is represented by f{circle around (x)}g, where:DFT[f{circle around (x)}g]=DFT[f]*DFT[g]=F(ω)G(ω).

The cross-correlation of two signals is represented by f*g, where:DFT[f*g]=DFT[ f]*DFT[g]= F(ω)G(ω).

The multiplication of two signals is represented by fg, where:DFT[fg]=DFT[f]{circle around (x)}DFT[g]=F(ω){circle around (x)}G(ω).

Projection theorem states:

$p_{\theta }\left(t\right)=\sum _m^M\sum _{n=1}^N[\hspace{1em}{h\left(m,n,t\right]}_{\theta }$ $P_{\theta }\left(\omega \right)=\mathrm{FT}\left[p_{\theta }\left(t\right)\right]=\mathrm{FT}\left\{\sum _m^M\sum _{n=1}^N{[h\left(m,n,t\right]\hspace{1em}}_{\theta }\right\}$ $P_{\theta }\left(\omega \right)=H\left(\omega \cos \theta ,\omega \sin \theta \right).$

In one embodiment of the present invention, the system is a patch of terrain, the input to this system is a chirp modulated on a microwave carrier frequency and the output of the system is the received backscattered 2D phase history data. Let h(m,n) be the impulse response, x(m,n) be the input (i.e., the emitted chirp), and y(m,n) (i.e., the received data) be the output of the system. Accordingly, y(m,n)=x(m,n){circle around (x)}h(m,n), or, in the frequency domain, Y(u,v)=X(u,v)H(u,v), where Y(u,v)=FFT[y(m,n)], X(u,v)=FFT[x(m,n)], and H(u,v)=FFT[h(m,n)].

The impulse response may be computed by

$H\left(u,v\right)=\frac{Y\left(u,v\right)}{X\left(u,v\right)}$ and then h(m,n)=FFT⁻¹[H(u,v)], noting that x(m,n), h(m,n), and y(m,n) are complex functions.

Also, the SAR system may be characterized by the following parameters: PRF (Pulse Repetition Frequency); RSF (Range Sampling Frequency); k (Chirp slope); τ_p (Pulse duration); λ (Radar wavelength); β (Real beam size); D (Real aperture size); β_s (Synthetic beam size); L_s (Synthetic aperture size); and v (Platform velocity).

To construct x(m,n) and form y(m,n) from the range time domain and azimuth spatial domain, first the pulse/chirp generated by the pulse generator is mathematically modeled. The generated pulse/chirp is a sinusoidal pulse with a Pulse Repetition Frequency equal to PRF, and may be represented as:

$x\left(t\right)=\left\{\Pi \left(\frac{t}{{\tau }_a}\right)*\exp \left(-i{\mathrm{kt}}^2\right)\right\},t<\frac{{\tau }_a}{2},$ where

$\Pi \left(\frac{t}{{\tau }_a}\right)$ is the two-dimensional RECT function:

$\Pi \left(\frac{t}{{\tau }_a}\right)=\{\begin{array}{c}1,t<\frac{{\tau }_a}{2}\\ 0.\end{array}$

Then the modulated chirp may be modeled as:

$x\left(t\right)\exp \left(-i{\omega }_{c}t\right)=\Pi \left(\frac{t}{{\tau }_p}\right)\exp \left[-\left(i{\omega }_{c}t+{\mathrm{kt}}^2\right)\right],t<\frac{{\tau }_p}{2}$ $\mathrm{FT}\left[x\left(t\right)\exp \left(-{\mathrm{i\omega }}_{c}t\right)\right]=X\left(\omega -{\omega }_c\right)$

This modulated chirp is emitted/transmitted by the antenna and subsequently illuminates a target in a region of interest or a terrain patch.

The number of pulses per second, or Pulse Repetition Frequency, is equal to PRF. The number of the samples in range direction (at a given θ), NR, is equal to the Listening Range Echo Time (LRET) multiplied by the Range Sampling Frequency:NR=LRET*RSF.

The number of the samples in the azimuth direction, NA, is equal to Listening Azimuth Echo Time (LAET) multiplied by the Pulse Repetition Frequency:

$\mathrm{NA}=\mathrm{LAET}*\mathrm{PRF}=\frac{L_s}{\nu }{\mathrm{PRF}}_t$ where v is the sensor velocity and L_s the systematic aperture length.

All the points in the field of the view are illuminated simultaneously while the chirp beam footprint (i.e., the beam spot size) interacts with the target. The terrain surface cover modifies the chirp (e.g., the amplitude, phase, and frequency), whereby the distorted and delayed chirp echoes back toward the antenna (i.e., as a backscattering signal).

The backscattering signal of a target point has an amplitude (i.e., signal strength) and a phase history. Since the SAR antenna is traveling, the echo of each point is received many times, or, in other words, the phase history of each target point is the phase delay of the point, which affects the phase delay of every other point in the field of view. Accordingly, the phase history of a point (k,l) in an M by N array of patch points for a given scene may be given by:y(k,l)=Ae^−i(nk+ml),0≦n≦N,0≦n≦N,0≦m≦M.

Each point on the ground is imaged many times as an airborne or spaceborne SAR platform moves along its path while measuring many points on the ground simultaneously.

To clarify, the received signal for an incidence angle/pulse may be modeled by:

$p_{\theta }\left(t\right)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}[\hspace{1em}\hspace{1em}\hspace{1em}{h\left(m,n,t\right]}_{\theta }=\sum _{m=0}^{N-1}\sum _{n=0}^{N-1}{\left(A_{\left(m,n\right)}{e}^{{\mathrm{i\varphi }}_{\left(m,n\right)}\left(t\right)}\right)}_{\theta },$ where A_(m,n) is the surface reflectance of the (m,n)_th target point and φ_(m,n)(t) is the phase history delay of the (m,n)_th point. The phase history delay is related to the time of an echo return

$\left(i.e.,\frac{2R_{\theta }}{c}\right).$

The received data is the convolution of the transmitted pulsed beams and target topographical points, where for one look, for example, a terrain patch is illuminated from many different incidence angles and is, thus, unfocused and must be de-convolved improve its resolution in both cross-track and along-track directions.

The data processing is first done on the range (i.e., cross-track) and then on the azimuth (i.e., along-track) directions.

In general, the received signal is given by:y(t)=x(t)exp(iω_ct){circle around (x)}p_θ(t)y(t)=FT⁻¹[(X(ω−ω_c)(P_θ(ω))]

For demodulation, a low pass filter may be used:

$\mathrm{LOW\_PASS}\mathrm{\_FILTER}=\Pi \left(\frac{\omega }{{\omega }_c}\right)\{\begin{array}{c}1,\omega <\frac{{\omega }_c}{2}\\ 0\end{array}$ such that

$X\left(\omega \right)=\left\{\Pi \left(\frac{\omega }{{\omega }_c}\right)*\left[X\left(\omega -{\omega }_c\right)\right]\right\},$ whereY(ω)=X(ω)P_θ(ω)=X(ω)H(ω cos θ, ω sin θ).

Tables 1–4, below, summarize the mathematical modeling, described above, of the aircraft position and beam incident angle, the transmitted signal, the preprocessed received signal in the time domain, and the Fourier transform of the preprocessed received signal for times t₀, t₁, t₂, and t_NA−1, respectively.

TABLE 1
Aircraft Position and beam   t0
incident angle at time t0    θ0 = θ(t0)
                             AP(t0)
Transmitted Signal           x(t0)exp(−iωct0) =
                             Π                      ⁡                                              (                                                                              t                            0                                                                                τ                            p                                                                          )                                                              ⁢                                          exp                      [                                              -                                                  (                                                                                    i                              ⁢                                                                                                                          ⁢                                                              ω                                c                                                            ⁢                                                              t                                0                                                                                      +                                                          kt                              0                              2                                                                                )                                                                    ]
Preprocessed Received Signal y(t0) = x(t0)   p(t0)
(time domain)                p                      ⁡                                              (                                                  t                          0                                                )                                                              =                                                                  ∑                                                  m                          =                          0                                                                          M                          -                          1                                                                    ⁢                                                                                          ⁢                                                                        ∑                                                      n                            =                            0                                                                                N                            -                            1                                                                          ⁢                                                                                                  ⁢                                                  [                                                      h                            (                                                          m                              ,                              n                              ,                                                              t                                0                                                                                      ]
                             y(t0) = {y(t0,0), y(t0,1), . . . , y(t0,NR−1)}
                             x(t0) = {x(t0,0), x(t0,1), . . . , x(t0,NR−1)}
                             p(t0) = {p(t0,0), p(t0,1), . . . , p(t0,NR−1)}
                             Δ                      ⁢                                                                                          ⁢                                              t                        r                                                              =                                                                                            t                                                      0                            ,                                                          i                              +                              1                                                                                                      -                                                  t                                                      0                            ,                            i                                                                                              =                                              1                        RSF
Fourier Transform            Y(ω0)
                             X(ω0)
                             P(ω0)
                             H(ω0cos(θ0), ω0sin(θ0))

TABLE 2
Aircraft Position and beam   t1
incident angle at time t1    Δta = t1 − t0
                             Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                                                  Δ                        ⁢                                                                                                  ⁢                        x                                            v                                                                                                                                                              Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                          1                      PRF
                             θ1 = θ(t1)
                             AP(t1)
Transmitted Signal           x(t1)exp(−iωct1) =
                             Π                      ⁡                                              (                                                                              t                            1                                                                                τ                            p                                                                          )                                                              ⁢                                          exp                      ⁡                                              [                                                  -                                                      (                                                                                          i                                ⁢                                                                                                                                  ⁢                                                                  ω                                  c                                                                ⁢                                                                  t                                  1                                                                                            +                                                              kt                                1                                2                                                                                      )                                                                          ]
Preprocessed Received Signal y(t1) = x(t1)   p(t1)
(time domain)                p                      ⁡                                              (                                                  t                          1                                                )                                                              =                                                                  ∑                                                  m                          =                          0                                                                          M                          -                          1                                                                    ⁢                                                                                          ⁢                                                                        ∑                                                      n                            =                            0                                                                                N                            -                            1                                                                          ⁢                                                                                                  ⁢                                                  [                                                      h                            (                                                          m                              ,                              n                              ,                                                              t                                1                                                                                      ]
                             y(t1) = {y(t1,0), y(t1,1), . . . , y(t1,NR−1)}
                             x(t1) = {x(t1,0), x(t0,1), . . . , x(t1,NR−1)}
                             p(t1) = {p(t1,0), p(t1,1), . . . , p(t1,NR−1)}
                             Δ                      ⁢                                                                                          ⁢                                              t                        r                                                              =                                                                                            t                                                      1                            ,                                                          i                              +                              1                                                                                                      -                                                  t                                                      1                            ,                            i                                                                                              =                                              1                        RSF
Fourier Transform            Y(ω1)
                             X(ω1)
                             P(ω1)
                             H(ω1cos(θ1), ω1sin (θ1))

TABLE 3
Aircraft Position and beam   t2
incident angle at time t2    Δta = t2 − t1
                             Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                                                  Δ                        ⁢                                                                                                  ⁢                        x                                            v                                                                                                                                                              Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                          1                      PRF
                             θ2 = θ(t2)
                             AP(t2)
Transmitted Signal           x(t2)exp(−iωct2) =
                             Π                      ⁡                                              (                                                                              t                            2                                                                                τ                            p                                                                          )                                                              ⁢                                          exp                      ⁡                                              [                                                  -                                                      (                                                                                          i                                ⁢                                                                                                                                  ⁢                                                                  ω                                  c                                                                ⁢                                                                  t                                  2                                                                                            +                                                              kt                                21                                2                                                                                      )                                                                          ]
Preprocessed Received Signal y(t2) = x(t2)   p(t2)
(time domain)                p                      ⁡                                              (                                                  t                          2                                                )                                                              =                                                                  ∑                                                  m                          =                          0                                                                          M                          -                          1                                                                    ⁢                                                                                          ⁢                                                                        ∑                                                      n                            =                            0                                                                                N                            -                            1                                                                          ⁢                                                                                                  ⁢                                                  [                                                      h                            (                                                          m                              ,                              n                              ,                                                              t                                2                                                                                      ]
                             y(t2) = {y(t2,0), y(t2,1), . . . , y(t2,NR−1)}
                             x(t2) = {x(t2,0), x(t2,1), . . . , x(t2,NR−1)}
                             p(t2) = {p(t2,0), p(t2,1), . . . , p(t2,NR−1)}
                             Δ                      ⁢                                                                                          ⁢                                              t                        r                                                              =                                                                                            t                                                      2                            ,                                                          i                              +                              1                                                                                                      -                                                  t                                                      2                            ,                            i                                                                                              =                                              1                        RSF
Fourier Transform            Y(ω2)
                             X(ω2)
                             P(ω2)
                             H(ω2cos(θ2), ω2sin(θ2))

TABLE 4
Aircraft Position and beam   tNA−1
incident angle at time tNA−1 Δta = tNA−1 − tNA−2
                             Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                                                  Δ                        ⁢                                                                                                  ⁢                        x                                            v                                                                                                                                                              Δ                      ⁢                                                                                          ⁢                                              t                        a                                                              =                                          1                      PRF
                             NA = Ls * FRF
                             θNA−1 = θ(tNA−1)
                             P(tNA−1)
Transmitted Signal           x(tNA−1)exp(−iωctNA−1) =
                             Π                      ⁡                                              (                                                                              t                                                          NA                              -                              1                                                                                                            τ                            p                                                                          )                                                              ⁢                                          exp                      ⁡                                              [                                                  -                                                      (                                                                                          i                                ⁢                                                                                                                                  ⁢                                                                  ω                                  c                                                                ⁢                                                                  t                                                                      NA                                    -                                    1                                                                                                                              +                                                              kt                                                                  NA                                  -                                  1                                                                2                                                                                      )                                                                          ]
Preprocessed Received Signal y(tNA−1) = x(tNA−1)   p(tNA−1)
(time domain)                p                      ⁡                                              (                                                  t                                                      NA                            -                            1                                                                          )                                                              =                                                                  ∑                                                  m                          =                          0                                                                          M                          -                          1                                                                    ⁢                                                                                          ⁢                                                                        ∑                                                      n                            =                            0                                                                                N                            -                            1                                                                          ⁢                                                                                                  ⁢                                                  [                                                      h                            (                                                          m                              ,                              n                              ,                                                              t                                                                  NA                                  -                                  1                                                                                                                      ]
                             y(tNA−1) = {y(tNA−1,0), y(tNA−1,1), . . . ,
                             y(t2,NA−1,NR−1)}
                             x(tNA−1) = {x(tNA−1,0), x(tNA−1,1), . . . ,
                             x(t2,NA−1,NR−1)}
                             p(tNA−1) = {p(tNA−1,0), p(tNA−1,1), . . . ,
                             p(t2,NA−1,NR−1)}
                             Δ                      ⁢                                                                                          ⁢                                              t                        r                                                              =                                                                                            t                                                                                    NA                              -                              1                                                        ,                                                          i                              +                              1                                                                                                      -                                                  t                                                                                    NA                              -                              1                                                        ,                            i                                                                                              =                                              1                        RSF
Fourier Transform            Y(ωNA−1)
                             X(ωNA−1)
                             P(ωNA−1)
                             H(ωNA−1cos(θNA−1), ωNA−1sin(θNA−1))

Next, the NA polar impulse frequency responses of the target may be computed as:

$\left\{\begin{array}{c}H\left({\omega }_0\cos \left({\theta }_0\right),{\omega }_0\sin \left({\theta }_0\right)\right)\\ H\left({\omega }_1\cos \left({\theta }_1\right),{\omega }_1\sin \left({\theta }_1\right)\right)\\ H\left({\omega }_2\cos \left({\theta }_2\right),{\omega }_2\sin \left({\theta }_2\right)\right)\\ \dots \\ \dots \\ H\left({\omega }_{\mathrm{NA}-1}\cos \left({\theta }_{\mathrm{NA}-1}\right),{\omega }_{\mathrm{NA}-1}\sin \left({\theta }_{\mathrm{NA}-1}\right)\right)\end{array}\right\},$ where for each impulse frequency response, (H(ω cos(θ), ω sin(θ)), there are NR polar sample points (frequency domain). Therefore, there are NA×NR sample points of the target polar impulse response in the frequency domain.

Now, 2D Cartesian impulse frequency response of the target impulse response may represented by H(u,v)=H(ω_k,ω_k), or, simply, H(k,l), where 0≦k≦NA−1 and 0≦l≦NA−1, and may be interpolated from the 2D polar impulse frequency response points:

$\left\{\begin{array}{c}H\left(0,1\right),H\left(0,2\right),\dots H\left(0,\mathrm{NR}-1\right)\\ H\left(1,1\right),H\left(1,2\right),\dots H\left(1,\mathrm{NR}-1\right)\\ \dots \\ \dots \\ H\left(\mathrm{NA},1\right),H\left(\mathrm{NA}-1,2\right),\dots H\left(\mathrm{NA}-1,\mathrm{NR}-1\right)\end{array}\right\}=\mathrm{INTERPOLATION}\left\{\begin{array}{c}H\left({\omega }_0\cos \left({\theta }_0\right),{\omega }_0\sin \left({\theta }_0\right)\right)\\ H\left({\omega }_1\cos \left({\theta }_1\right),{\omega }_1\sin \left({\theta }_1\right)\right)\\ H\left({\omega }_2\cos \left({\theta }_2\right),{\omega }_2\sin \left({\theta }_2\right)\right)\\ \dots \\ \dots \\ H\left({\omega }_{\mathrm{NA}-1}\cos \left({\theta }_{\mathrm{NA}-1}\right),{\omega }_{\mathrm{NA}-1}\sin \left({\theta }_{\mathrm{NA}-1}\right)\right)\end{array}\right\}$

Finally the SAR complex image may be reconstructed by:

$h\left(m,n\right)={\mathrm{FFT}}^{-1}\left[H\left(u,v\right)\right]$ $h\left(m,n\right)=A\left(m,n\right)e^{i\varphi \left(m,n\right)}=h_{\mathrm{real}}\left(m,n\right)+{\mathrm{ih}}_{\mathrm{img}}\left(m,n\right)$ $A\left(m,n\right)=h\left(m,n\right)=\sqrt{h_{\mathrm{real}}^2\left(m,n\right)+h_{\mathrm{img}}^2\left(m,n\right)}$ $\varphi \left(m,n\right)=\arctan \frac{h_{\mathrm{img}}\left(m,n\right)}{h_{\mathrm{real}}\left(m,n\right)}$

Accordingly, the SAR phase history data, H(u,v), may be computed from the reconstructed SAR image, h(m, n), by taking the Fourier transform: H(u,v)=FFT [h(m,n)]. In one embodiment of the invention, only the magnitude of h(m,n) may be retained for display as the SAR image.

Referring now to the drawing, in which like reference numbers refer to like elements throughout the various figures that comprise the drawing, FIG. 3 illustrates plot 32 of complex phase history data, where plot 32 includes phase history data, H(u,v), plotted on the slow time (i.e., azimuth, along-track, cross-range) and fast time (i.e., range, cross-track) axes. The entire set of phase history data may be passed to inverse fast Fourier transformer (IFFT) 35 to form high resolution SAR image 37. Alternately, only phase history data in sub-aperture 34 may be passed to IFFT 36 to form low resolution SAR image 38.

Mathematically, taking a sub-aperture of the phase history data corresponds to multiplying the entire phase history data set, H(u,v), by a two-dimensional filter. For a rectangular sub-aperture, for example, this corresponds to multiplying by a 2D RECT function; for a circular sub-aperture, this corresponds to multiplying by a 2D Gaussian distribution.

A 1D box filter, for example, may be represented as:

$h\left(x\right)=\mathrm{RECT}\left(\frac{x}{\alpha }\right)=\Pi \left(\frac{x}{\alpha }\right)=\{\begin{array}{c}1,x<\frac{\alpha }{2}\\ 0\end{array}$ where taking the Fourier transform gives:FT(h(x))=F(u)=α sinc(αu).

A 2D box filter (i.e., a rectangle) may be represented as:

$g\left(x,y\right)=\mathrm{RECT}\left(\frac{x}{\alpha }\right)\mathrm{RECT}\left(\frac{y}{\beta }\right),$ where taking the Fourier transform gives:FT(g(x,y))=F(u,v)=α sinc(αu)β sinc(βv).

FIG. 4A illustrates a number of sub-apertures, 41–43, that may be taken on the phase history data set 40. Since each data point in the phase history data set 40 is characterized by a slow-time value, a fast-time value, and a magnitude and phase, multiplying the phase history data set 40 by any one of sub-apertures 41–43 corresponds to multiplying by a 2D RECT sub-aperture function, described above, which is a 2D slice of the 3D RECT sub-aperture 44 shown in FIG. 4B.

A 1D Gaussian filter, on the other hand, may be represented as:

$h\left(x\right)=\frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{{\left(x\right)}^2}{2{\sigma }^2}}$ where taking the Fourier transform gives:

$\mathrm{FT}\left[h\left(x\right)\right]=F\left(u\right)=e^{-\frac{{{\sigma }^2\left(u\right)}^2}{2}}.$

A 2D circularly symmetric Gaussian distribution may be represented as:

$h\left(x,y\right)=\frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{{\left(x\right)}^2}{2{\sigma }^2}}\frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{{\left(y\right)}^2}{2{\sigma }^2}}$ where the Fourier transform gives:

$F\left(u,v\right)=e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}.$

A 2D circularly asymmetric Gaussian distribution may be represented as:

$h\left(x,y\right)=\frac{1}{\sqrt{2\pi {\sigma }_x}}{e}^{-\frac{{\left(x\right)}^2}{2{\sigma }_x^2}}\frac{1}{\sqrt{2\pi {\sigma }_y}}{e}^{-\frac{{\left(y\right)}^2}{2{\sigma }_y^2}}$ where the Fourier transform gives:

$F\left(u,v\right)=e^{-\frac{{\sigma }_x^2{u}^2+{\sigma }_y^2{v}^2}{2}}.$

FIG. 5A illustrates a number of sub-apertures, 51–53, that may be taken on the phase history data set 50. Since each data point in the phase history data set 50 is characterized by a slow-time value, a fast-time value, and a magnitude and phase, multiplying the phase history data set 50 by any one of sub-apertures 51–53 corresponds to multiplying by a 2D circularly symmetric Gaussian sub-aperture function, described above, which is a 2D slice of the 3D Gaussian sub-aperture 55 shown in FIG. 5B.

According to other embodiments of the invention, various other sub-aperture geometries may be used in addition to the rectangular, circularly symmetric, and circularly asymmetric embodiments described above. FIG. 6 illustrates one such exemplary embodiment, where sectoral sub-apertures 61–64, for example, are taken on phase history data set 60.

Those skilled in the art will recognize the following properties of the Fourier Transform with respect to a 1D convolution:g(x)=f(x){circle around (x)}h(x)FT(g(x))=FT(f(x))FT(h(x))and with respect to a 2D convolution:g(x,y)=f(x,y){circle around (x)}h(x,y)FT(g(x,y))=FT(f(x,y))FT(h(x,y)

BY representing a high-quality reconstructed zero-frequency center SAR phase history data as F(u,v), a SAR image may be formed by taking the inverse Fourier transform (f(x,y)=FT⁻¹[F(u,v)]) or the inverse Fast Fourier transform (IFFT) (f(x,y)=FFT⁻¹[F(u,v)]).

In one embodiment, phase history data F(u,v) is filtered by multiplication with a 2D RECT function representing a rectangular sub-aperture, which corresponds to a convolution of the SAR image with the Fourier transform of the 2D RECT function, as described above. This is shown by:

$F\left(u,v\right)\mathrm{RECT}\left(\frac{u}{\alpha }\right)\mathrm{RECT}\left(\frac{v}{\beta }\right)f\left(x,y\right)\otimes \left[\mathrm{\alpha sin}c\left(\alpha x\right)\mathrm{\beta sin}c\left(\beta y\right)\right],$ where α and β determine the size of a rectangle in the Fourier domain (i.e., the phase history domain).

In another embodiment, phase history data F(u,v) is filtered by multiplication with a 2D circularly symmetric Gaussian function representing a circular sub-aperture, which corresponds to a convolution of the SAR image with the Fourier transform of the 2D circularly symmetric Gaussian function, as described above. This is shown by:

$F\left(u,v\right)e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}f\left(x,y\right)\otimes \left[\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x\right)}^2}{2{\sigma }^2}}\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(y\right)}^2}{2{\sigma }^2}}\right],$ where σ determines the size of the circle in the Fourier domain (i.e., the phase history domain).

In yet another embodiment, phase history data F(u,v) is filtered by multiplication with a 2D circularly asymmetric Gaussian function representing an elliptical sub-aperture, which corresponds to a convolution of the SAR image with the Fourier transform of the 2D circularly asymmetric Gaussian function, as described above. This is shown by:

$F\left(u,v\right)e^{-\frac{{\sigma }_x^2{u}^2+{\sigma }_y^2{v}^2}{2}}f\left(x,y\right)\otimes \left[\frac{1}{\sqrt{2\pi }{\sigma }_x}{e}^{-\frac{{\left(x\right)}^2}{2{\sigma }_x^2}}\frac{1}{\sqrt{2\pi }{\sigma }_y}{e}^{-\frac{{\left(y\right)}^2}{2{\sigma }_y^2}}\right],$ where σ_x and σ_y determine the size of the ellipse in the Fourier Domain (i.e., the phase history domain).

Accordingly, an embodiment of the present invention may iteratively filter the phase history data by iteratively multiplying the phase history data, in the Fourier domain, with 2D filter functions corresponding to different sub-apertures. Furthermore, the inverse Fourier transform of the filtered phase history data may be taken after each iteration in order to obtain a multi-resolution SAR image.

For example, a first sub-aperture filtered phase history data set may correspond to

$F\left(u,v\right)e^{-\frac{{\sigma }_1^2\left(u^2+v^2\right)}{2}},$ where the first SAR image may be obtained by taking the inverse Fourier transform:

${\mathrm{FT}}^{-1}\left\{F\left(u,v\right)e^{-\frac{{\sigma }_1^2\left(u^2+v^2\right)}{2}}\right\}.$ In some instances, the first SAR image may not be satisfactory for a number of reasons, which may include, for example, undesirable aliasing effects in the first SAR image that arise as a result of the first sub-aperture selected. Accordingly, a new sub-aperture may be chosen, where the new sub-aperture filtered phase history data set may correspond to

$F\left(u,v\right)e^{-\frac{{\sigma }_2^2\left(u^2+v^2\right)}{2}}.$ The new SAR image, therefore, may be obtained by taking the inverse Fourier transform:

$F^{-1}\left\{F\left(u,v\right)e^{-\frac{{\sigma }_2^2\left(u^2+v^2\right)}{2}}\right\}.$ This iterative process may be repeated by changing σ until a desirable SAR image is formed from the sub-aperture filtered phase history data set, and/or until a correlated image pyramid is formed for image compression, to be described below. In general:

$\left[F_{n+1}\left(u,v\right)\right]=\left\{\left[F_n\left(u,v\right)\right]e^{-\frac{{\sigma }_n^2\left(u^2+v^2\right)}{2}}\right\},\left[f_n\left(x,y\right)\right]={\mathrm{FT}}^{-1}\left\{F_n\left(u,v\right)\right\}$ where n is the iteration number.

In an alternate embodiment, the size of the filter is fixed, but the filter may be applied successively. For example, in the first iteration, the phase history data is filtered once:

$F\left(u,v\right)e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}},$ giving SAR image:

${\mathrm{FT}}^{-1}\left\{F\left(u,v\right)e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}\right\}.$

In the next iteration, the same filter is applied to the already once-filtered phase history data:

$\left\{F\left(u,v\right)e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}\right\}{e}^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}},$ giving SAR image:

${\mathrm{FT}}^{-1}\left\{F\left(u,v\right)e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}{e}^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}\right\}.$

In general:

$\left[F_{n+1}\left(u,v\right)\right]=\left\{\left[F_n\left(u,v\right)\right]e^{-\frac{{\sigma }^2\left(u^2+v^2\right)}{2}}\right\},\left[f_n\left(x,y\right)\right]={\mathrm{FT}}^{-1}\left\{F_n\left(u,v\right)\right\}$ where n is the iteration number.

Those skilled in the art will recognize that the iteratively generated SAR images, obtained by exemplary embodiments of the invention such as the ones described above, may be spatially correlated in the image domain. In one embodiment of the invention, therefore, obtaining a plurality of correlated SAR images from different sub-apertures may include generating a correlated image pyramid. Accordingly, a Laplacian pyramid technique may be used to compress the plurality of correlated SAR images in the correlated image pyramid. In a further embodiment, target detection and surface or subsurface type classification may be performed by taking the integral of each level in the correlated image pyramid, where each level may correspond to a distinct sub-aperture filtered phase history data set.

In another embodiment, the phase history data may also be iteratively smoothed by: F_n+1(u,v)={F_n(u,v){circle around (x)}H(u,v)}, where H(u,v) is the Fourier transform of an optimal filter and n is the iteration number. As described above, convolution in one domain corresponds to multiplication in the transform domain, thus f_n+1(x,y)=f_n(x,y)h(x, y).

In one embodiment of the invention, a plurality of correlated SAR images may be obtained by iteratively changing a sub-aperture on the phase history data, as described above. In another embodiment of the invention, a plurality of correlated SAR images may be obtained by successively applying a sub-aperture on previously filtered phase history data, as described above. In yet another embodiment, a plurality of correlated SAR images may be obtained by choosing one particular SAR image and then successively downsampling it (e.g., successively convolving the image with a weighting function). Accordingly, the plurality of SAR images that are obtained may be compressed using a Laplacian pyramid, for example, by taking advantage of the correlation between each image in the plurality of images.

In an embodiment of the invention, therefore, successively downsampling a SAR image to form a plurality of SAR images may begin by setting the frame as the entire phase history data set, ph₀=F(u,v), or, alternately, by choosing a sub-aperture on the phase history data set, ph₀=Sub_Aperture[F(u,v)], where F(u,v) represents the phase history data and choosing a sub-aperture on the phase history data set corresponds to multiplying the phase history data set by a 2D function representing the sub-aperture, as described above. Then, to obtain a lower resolution SAR image for each i^th iteration, the higher resolution SAR image from the previous (i−1)^th iteration is low-pass filtered to produce the lower resolution SAR image (i.e., F^T−1(ph_i)=LPF[FT⁻¹(ph_i)]). This corresponds to convolving the image with a weighting function (such as a Gaussian weighting function, for example) in the spatial domain, or, equivalently, multiplying the image with the weighting function in the frequency domain. Accordingly, a sequence of multi-resolution images {FT⁻¹(ph₀), FT⁻¹(ph₁), . . . , FT⁻¹(ph_n) } are generated, where the image generated by each i^th iteration may be a fourth of the size of the image of the previous (i−1)^th iteration

$\left(i.e.,\mathrm{Size\_of}\left({\mathrm{ph}}_i\right)=\frac{1}{4}\mathrm{Size\_of}\left({\mathrm{ph}}_{i-1}\right)\right),$ depending on the weighting function chosen. Those skilled in the art will recognize that a plurality of different weighting functions such as a Gaussian, triangular, rectangular, trimodal, broad Gaussian, or other weighting function may be chosen without departing from the scope of the present invention.

FIG. 7A is a pictorial illustrating the effect of successively downsampling (i.e., low-pass filtering) a 1D complex image (g₀) having complex data points 701–709 in order to obtain lower resolution complex images g₁ and g₂ in a correlated image pyramid, according to an exemplary embodiment of the present invention. The exemplary embodiment utilizes the exemplary Gaussian weighting functions illustrated in FIGS. 7B and 7C. Gaussian weighting function 73, therefore, is centered on point 703 in image g₀, and convolved with image g₀ in order to obtain a value for point 710 in image g₁ that corresponds to a weighted sum of the values of points 701–705. Gaussian function 75 is centered on point 705 in image g₀ and is convolved with image g₀ in order to obtain a value for point 711 in image g₂ that corresponds to a weighted sum of the values of points 703–707. Gaussian function 77 is centered on point 707 and is convolved with image g₀ in order to obtain a value for point 712 in image g₁ that corresponds to a weighted sum of the values of points 705–709. To obtain the next lower resolution image, g₂, Gaussian function 71 is centered on point 712 in image g₁, and convolved with image g₁ in order to obtain a value for point 720 in image g₂ that corresponds to a weighted sum of the values of points 710–714. Images g₀, g₁, and g₂ are, therefore, correlated. Accordingly, g₁ may be subtracted from g₀ (giving L₀) and g₂ may be subtracted from g₁ (giving L₁) in order to obtain a Laplacian pyramid (i.e., L₀, L₁, and g₂ as L₂), which may encoded to provide a compressed representation of 1D images g₀, g₁, and g₂.

In general, by downsampling (i.e., low-pass filtering) the SAR images successively with a low-pass filter, high frequency content (i.e., edges) is removed from the image. As such, a Laplacian pyramid of edge enhanced images may be generated by subtracting the SAR image of each (i+1)^th iteration from the SAR image of the previous i^th iteration (i.e., L_i=FT⁻¹(ph_i)−FT⁻¹(ph_i+1)). Accordingly, a sequence of Laplacian edge-enhanced images, L₀, L₁, . . . , L_n, may be generated and subsequently encoded to achieve greater data compression, where each L_i is smaller than its predecessor, L_i−1, by a scale of one half due to reduced image sample density in the horizontal and vertical axes.

FIG. 8 is a block diagram illustrating an exemplary method of obtaining a compressed representation of multi-resolution SAR images g₀–g₃, according to one embodiment of the present invention. In such an embodiment, it is assumed that multi-resolution SAR images g₀–g₃ are correlated, thereby constituting a correlated image pyramid. Accordingly, high resolution SAR image g₀ may be subtracted by first lower resolution SAR image g₁ by subtractor 803 to obtain first Laplacian image L₀; first lower resolution SAR image g₁ may be subtracted by second lower resolution SAR image g₂ by subtractor 802 to obtain second Laplacian image L₁; and second lower resolution SAR image g₂ may be subtracted by lowest resolution SAR image g₃ by subtractor 801 to obtain third Laplacian image L₂; and lowest resolution SAR image g₃ may be used as fourth Laplacian image L₃. Each of Laplacian images L₀–L₃, constituting a Laplacian image pyramid, may then be encoded as quantized values Q₀–Q₃, respectively. Quantized values Q₀–Q₃, therefore, provide a compressed representation of multi-resolution SAR images g₀–g₃ in the correlated image pyramid. Those skilled in the art will recognize that a greater correlation between respective SAR image levels provides greater compression.

Reconstruction of original SAR images g₀–g₃ may begin by obtaining quantized value Q₃, which is then decoded to obtain a reconstructed fourth Laplacian image L₃* (which is also equivalent to reconstructed lowest resolution SAR image g₃*). Quantized value Q₂ is obtained next and subsequently decoded to obtain reconstructed third Laplacian image L₂*. Reconstructed third Laplacian image L₂* may then be added to reconstructed lowest resolution SAR image g₃ by adder 804 to obtain reconstructed second lower resolution SAR image g₂*. Quantized value Q₁ is obtained next and subsequently decoded to obtain reconstructed second Laplacian image L₁. Reconstructed second Laplacian image L₁ may then be added to reconstructed second lower resolution SAR image g₂* by adder 805 to obtain reconstructed first lower resolution SAR image g₁. Finally, quantized value Q₀ is obtained and subsequently decoded to obtain reconstructed first Laplacian image L₀*. Reconstructed first Laplacian image L₀* may then be added to reconstructed first lower resolution SAR image 91 by adder 806 to obtain reconstructed high resolution SAR image g₀*. Those skilled in the art will recognize the reconstruction steps described above may constitute a progressive transmission of multi-resolution SAR images to a client. Accordingly, such a progressive transmission may stop whenever the client has received a SAR image having an adequate resolution, thereby possibly precluding the need to transmit all of quantized values Q₀–Q₃.

Those skilled in the art will recognize that many types of encoding schemes may be used to obtain quantized values without departing from the scope of the present invention, including, for example, entropy encoding techniques such as Huffman coding, range encoding, and arithmetic encoding, or static coding techniques such as Unary coding, Elias coding, Fibonacci coding, Golomb coding, or Rice coding.

FIG. 9 is a flow diagram illustrating an exemplary method for forming a correlated image pyramid and a corresponding Laplacian pyramid for compression. In step 901, chirps/pulses modulated on a microwave sinusoidal carrier signal are transmitted from the SAR antenna of an airborne or spaceborne SAR imaging system. The chirp/pulses are subsequently reflected back toward the SAR imaging system from objects and terrain in the area being imaged. In step 902, the SAR antenna receives these modulated return signals (i.e., backscattered return pulses). Step 903 subsequently reconstructs a high quality phase history data set from the modulated return signals, as described above, and sets the iteration number, n, to “0”. High resolution SAR image 905 (i.e., G₀) may then be formed in step 904 by taking the IFFT of the entire phase history data set. Accordingly, in one embodiment (not shown in FIG. 9), high resolution SAR image 905 may be used to obtain a correlated image pyramid by successively downsampling (i.e., low-pass filtering) high resolution SAR image 905 and subtracting successive levels of the image pyramid to form a corresponding Laplacian pyramid, as described above. Alternately, step 906 may proceed to choose a first sub-aperture (SA₀) (e.g., rectangle, circle, ellipse, etc.) on the phase history data set, where the first sub-aperture is represented by a 2D function, f(SA₀). Step 907 then obtains the next image in the correlated image pyramid by taking the IFFT of the sub-aperture data set (i.e., by taking the IFFT of the product of the phase history data set and the 2D sub-aperture function, f(SA₀)), and then increments the iteration number, n, to “1”. The remainder of the correlated image pyramid may be obtained according to a plurality of embodiments of the present invention.

In an embodiment of the invention, the method of forming the correlated image pyramid, described above, may be iteratively repeated if any images in the pyramid have undesirable characteristics. Such undesirable characteristics may include, but are not limited to, image distortion, aliasing, blur, noise, and the like.

In one embodiment, for example, the path labeled “Option A” in FIG. 9 may be taken, where step 910 obtains a first Laplacian image, L₀, by subtracting the SAR image of the current iteration, G₁, from the SAR image of the previous iteration, G₀. In step 911, the SAR image of the current iteration, G₁, is downsampled (i.e., low-pass filtered), as previously described, to obtain the SAR image of the next iteration, G₂. Step 912 then determines whether the current iteration number, n, is less than or equal to a total number of desired iterations, N, and then increments the current iteration number, n. If n is less than or equal to N, the process moves back to step 910 to form the next Laplacian image, and then to step 911 to form the next SAR image. Accordingly, steps 910–912 are repeated until the correlated image pyramid and corresponding Laplacian pyramid are formed for a desired number of iterations (or levels), N. Once the desired number of iterations are completed, the process ends in step 900.

FIG. 10 illustrates the formation of a correlated image pyramid and corresponding Laplacian pyramid for four iterations (or levels), according to the embodiment of the invention described above with respect to “Option A”. Accordingly, entire phase history data set 102 is inverse Fourier transformed by IFFT 105 to obtain high resolution SAR image 107, which corresponds to highest resolution image g₀ in the correlated image pyramid. Sub-aperture filtered phase history data set 104 is also inverse Fourier transformed, by IFFT 106, to form first lower resolution SAR image 108, which corresponds to first lower resolution image g₁ in the correlated image pyramid. First lower resolution image g₁ may then be downsampled (i.e., low-pass filtered), as described above, to obtain second lower resolution image g₂. Second lower resolution image g₂ may be subsequently downsampled to form lowest resolution image g₃. The corresponding Laplacian pyramid may be formed by subtracting g₁ from g₀ in subtractor 1000 to form first Laplacian image L₀, g₂ from g₁ in subtractor 1001 to form second Laplacian image L₁, and g₃ from g₂ in subtractor 1002 to form third Laplacian image L₂. Lowest resolution image g₃ may be used as fourth Laplacian image L₃. Those skilled in the art will recognize that a sub-aperture of any desirable size and/or shape may be chosen, that the highest resolution image g₀ may be taken from any sub-aperture filtered phase history data set, and that any number of iterations may be performed to obtain a desired number of levels of the correlated image pyramid.

In another embodiment, the path labeled “Option B” in FIG. 9 may be taken, where step 920 obtains a first Laplacian image, L₀, by subtracting the SAR image of the current iteration, G₁, from the SAR image of the previous iteration, G₀. In step 921, the SAR image of the current iteration, G₁, is convolved with a low-pass filter corresponding to the first sub-aperture. In other words, the phase history data corresponding to the SAR image of the current iteration (i.e., FFT(G₁)), is multiplied by the 2D sub-aperture function, f(SA₀). Step 922 obtains the resulting SAR image of the next iteration, G₂, by taking the inverse Fourier transform of the twice filtered phase history data set:FT⁻¹{[PH·f(SA₀)]·f(SA₀)}=FT⁻¹{FT(G₁)·f(SA₀)}.

Step 923 then determines whether the current iteration number, n, is less than or equal to a total number of desired iterations, N, and then increments the current iteration number, n. If n is less than or equal to N, the process moves back to step 920 to form the next Laplacian image, and then to steps 921–922 to form the next SAR image. Accordingly, steps 920–922 are repeated until the correlated image pyramid and corresponding Laplacian pyramid are formed for a desired number of iterations (or levels), N. Once the desired number of iterations are completed, the process ends in step 900.

FIG. 11 illustrates the formation of a correlated image pyramid and corresponding Laplacian pyramid for four iterations (or levels), according to the embodiment of the invention described above with respect to “Option B”. Accordingly, phase history data set 1100 is inverse Fourier transformed by IFFT 1110 to obtain highest resolution SAR image g₀ in the correlated image pyramid. The first lower resolution SAR image g₁ in the correlated image pyramid may be obtained by filtering phase history data set 1100 by a circularly symmetric sub-aperture. Accordingly, sub-aperture data set 1103 is inverse Fourier transformed, by IFFT 1113, to form first lower resolution SAR image g₁. The second lower resolution SAR image g₂ in the correlated image pyramid may be obtained by filtering sub-aperture data set 1103 by the same circularly symmetric sub-aperture. Accordingly, sub-aperture data set 1102 is inverse Fourier transformed, by IFFT 1112, to form second lower resolution SAR image g₂. The lowest resolution SAR image g₃ in the correlated image pyramid may be obtained by filtering sub-aperture data set 1102 by the same circularly symmetric sub-aperture. Accordingly, sub-aperture filtered phase history data set 1101 is inverse Fourier transformed, by IFFT 1111, to form lowest resolution SAR image g₃.

The corresponding Laplacian pyramid may then be formed by subtracting g₁ from g₀ in subtractor 1105 to form first Laplacian image L₀, g₂ from g₁ in subtractor 1106 to form second Laplacian image L₁, and g₃ from g₂ in subtractor 1107 to form third Laplacian image L₂. Lowest resolution image g₃ may be used as fourth Laplacian image L₃. Those skilled in the art will recognize that a sub-aperture of any desirable size and/or shape may be chosen, that the highest resolution image g₀ may be taken from any sub-aperture filtered phase history data set, and that any number of iterations may be performed to obtain a desired number of levels of the correlated image pyramid.

In yet another exemplary embodiment, the path labeled “Option C” in FIG. 9 may be taken, where step 930 obtains a first Laplacian image, L₀, by subtracting the SAR image of the current iteration, G₁, from the SAR image of the previous iteration, G₀. In step 931, a new sub-aperture (SA₁) (e.g., rectangle, circle, ellipse, etc.) on the phase history data set is chosen for the current iteration, where the new sub-aperture is represented by a 2D function, f(SA₁). Step 932 then obtains the next image, G₂, in the correlated image pyramid by taking the IFFT of the sub-aperture data set (i.e., by taking the IFFT of the product of the phase history data set and the 2D sub-aperture function, f(SA₁)), and then increments the iteration number, n, to “2”. Step 933 then determines whether the current iteration number, n, is less than or equal to a total number of desired iterations, N, and then increments the current iteration number, n. If n is less than or equal to N, the process moves back to step 930 to form the next Laplacian image, and then to steps 931–932 to choose a new sub-aperture and form the next SAR image. Accordingly, steps 930–932 are repeated until the correlated image pyramid and corresponding Laplacian pyramid are formed for a desired number of iterations (or levels), N. Once the desired number of iterations are completed, the process ends in step 900.

FIG. 12 illustrates the formation of a correlated image pyramid and corresponding Laplacian pyramid for four iterations (or levels), according to the embodiment of the invention described above with respect to “Option C”. Accordingly, phase history data set 120 is inverse Fourier transformed by IFFT 1210 to obtain highest resolution SAR image g₀ in the correlated image pyramid. The first lower resolution SAR image g₁ in the correlated image pyramid may be obtained by filtering phase history data set 120 by a rectangular sub-aperture. Accordingly, sub-aperture data set 121 is inverse Fourier transformed, by IFFT 1211, to form first lower resolution SAR image g₁. The second lower resolution SAR image g₂ in the correlated image pyramid may be obtained by filtering phase history data set 120 by a new rectangular sub-aperture. Accordingly, sub-aperture data set 122 is inverse Fourier transformed, by IFFT 1212, to form second lower resolution SAR image g₂. The lowest resolution SAR image g₃ in the correlated image pyramid may be obtained by filtering phase history data set 120 by yet another new rectangular sub-aperture. Accordingly, sub-aperture data set 123 is inverse Fourier transformed, by IFFT 1213, to form lowest resolution SAR image g₃.

The corresponding Laplacian pyramid may then be formed by subtracting g₁ from g₀ in subtractor 1200 to form first Laplacian image L₀, g₂ from g₁ in subtractor 1201 to form second Laplacian image L₁, and g₃ from g₂ in subtractor 1202 to form third Laplacian image L₂. Lowest resolution image g₃ may be used as fourth Laplacian image L₃. Those skilled in the art will recognize that new sub-apertures of any desirable size and/or shape may be chosen, that the highest resolution image g₀ may be taken from any sub-aperture filtered phase history data set, and that any number of iterations may be performed to obtain a desired number of levels of the correlated image pyramid.

In a further embodiment of the present invention, the path labeled “Target Detection” in FIG. 9 may be taken before, after, or in parallel with one or more of the other paths (i.e., “Option A”, “Option B”, and “Option C”). In step 991, a first sub-aperture filtered data set is formed by multiplying the phase history data set with a function representing the current sub-aperture, f(SA₀). By integrating over the first sub-aperture filtered data set, step 991 obtains a first Feature, F₀, which is stored in step 992. Step 993 then determines whether the current iteration number, n, is less than or equal to a total number of desired iterations, N, and then increments the current iteration number, n. If n is less than or equal to N, the process moves to step 990 to choose a new sub-aperture, and then to steps 991–992 to form and store the feature, F_n, for the current iteration. Accordingly, steps 990–992 are repeated until the desired number of features, N, are formed and stored. Once the desired number of iterations are completed, the process moves to step 994 to detect a target based on these N stored features. Target detection may include using the N features as input to a supervised (Maximum Likelihood), a Neural Networks or unsupervised classifier such as a K-means or ISODATA clustering method for target, texture, and geometric cover type classifications. If, in step 995, the Receiver Operating Characteristic (ROC) for the current target detection is acceptable, then target detection is completed in step 999. If the ROC is unacceptable, the iteration number, n, is reset and the target detection process restarts from step 990. Accordingly, an embodiment of the present invention provides iterative target detection on a captured scene characterized by phase history data, by iteratively choosing new sub-apertures on the phase history data for forming features to be used in target detection algorithms.

FIG. 13 illustrates the formation of four features for use in target detection, according to the embodiment of the invention described above with respect to “Target Detection”. Accordingly, phase history data within first, second, third, and fourth sectoral sub-apertures 1301–1304 are integrated by integrators 1311–1314, respectively. The resulting four features are inputted to a statistical (maximum likelihood) or neural network (packpropagation) classifier 1300. Other inputs to classifier 1300 include prototype vectors for a plurality of targets or classes (i.e., prototype vector for target or class 1 to N).

Exemplary features may also include annuli and sector masks, where each annuli and sector corresponds to a sub-aperture. Accordingly, in a target detection method according to the present invention, only a few concentric annuli in an upper Fourier half-plane (SAR phase history data) and equal-angle sectors in a lower Fourier half-plane, for selected hyper spectral bands or in an orthogonal subspace, features may be extracted.

In yet another embodiment of the invention, multi-resolution SAR images may be progressively transmitted to a client that requests a desired image of a scene or one or more targets within the scene. The client may require that the desired image have a desired resolution and be substantially free from undesirable characteristics such as aliasing, noise, blur, and the like, for example. In general, the method may begin by obtaining SAR phase history data corresponding to the scene. A correlated image pyramid may then be generated, as described above, comprising a plurality of multi-resolution images of the scene. A Laplacian pyramid may then be generated, as described above, comprising a plurality of edge-enhanced images corresponding to a difference in the multi-resolution images in successive levels of the correlated image pyramid. Each one of the Laplacian pyramid images may then be encoded for further compression.

A client request for the desired image may be received either before, after, or simultaneous to the steps described above. Once the encoded Laplacian pyramid images are ready, they are successively transmitted to the client, starting from the bottom level of the pyramid (i.e., lowest resolution and smallest size image) and working up the pyramid until the desired image has been sent.

In one embodiment of the invention, a request for a desired image is received from a client, but the characteristics of the desired image are unknown at the time of transmission. Accordingly, the method may continue to transmit the encoded images starting from a bottom level until an indication is received from the client that the desired image has been transmitted. The method may stop sending images once it's received such an indication. In a further embodiment, the method may transmit all the encoded images without receiving an indication that the desired image has been transmitted. This corresponds to a situation where the client is not happy with the characteristics of the images in the correlated image pyramid. In such a case, therefore, a new correlated image pyramid and corresponding encoded Laplacian pyramid may be formed, and the method repeated, until the desired image is transmitted.

In another embodiment, the method may also detect one or more targets in the scene, as described above. Accordingly, correlated image pyramids may be generated as comprising a plurality of multi-resolution images of one or more desired targets, rather than the entire scene. The method may then receive a request from a client for a desired image of a particular target in the scene. The appropriate image pyramid is then selected and the progressive transmissions commenced until the desired image is transmitted.

Although illustrated and described above with reference to certain specific embodiments, the present invention is nevertheless not intended to be limited to the details shown. Rather, various modifications may be made in the details within the scope and range of equivalents of the claims and without departing from the invention.

Claims

1. A method of generating multi-resolution synthetic aperture radar (SAR) images, the method comprising the steps of: a) obtaining SAR phase history data, denoted as PH;b) calculating an inverse Fourier transform of a first data set, denoted as ph0, of the SAR phase history data, PH, to form a highest resolution image, denoted as G0, where G0=FT−1(ph0);c) calculating the inverse Fourier transform of a plurality of other data sets, denoted as ph1 . . . phN, to form a plurality of correlated lower resolution images, denoted as G1 . . . GN, where Gi=FT−1(phi), for i=1 to N; d) subtracting successive ones of the highest resolution image, G0, and the lower resolution images, G1 . . . GN, to obtain a plurality of edge enhanced images, denoted as L0 . . . LN−1, where Li=Gi−Gi+1, for i=0 to N−1; ande) obtaining a lowest resolution image, denoted as LN.

2. A method according to claim 1, wherein steps (b)–(f) are repeated for a new first data set and a new plurality of other data sets, if one or more of the highest resolution image and the plurality of correlated lower resolution images have undesirable characteristics.

3. A method according to claim 1, further including the step of encoding the plurality of edge enhanced images and the lowest resolution image.

4. The method of claim 3, wherein the encoding step encodes according to one of a Huffman code, a range code, an arithmetic code, a Unary code, an Elias code, a Fibonacci code, a Golomb code, or a Rice code.

5. The method of claim 1, wherein the first data set, ph0, comprises the SAR phase history data, PH: ph0=PH.

6. The method of claim 1, wherein the first data set, ph0, comprises a first sub-aperture set, denoted as SA0, of the SAR phase history data, PH, and step (b) further comprises multiplying the SAR phase history data, PH, by a function representing the first sub-aperture set, denoted as fSA0, where ph0=PH·fSA0.

7. The method of claim 6, wherein the first sub-aperture set, SA0, is chosen from one of a circularly symmetric Gaussian, a circularly asymmetric Gaussian, a rectangle, an annuli, or a sector.

8. A method according to claim 6, further comprising the step of successively multiplying the first data set, ph0, by the function representing the first sub-aperture set, fSA0, to obtain the plurality of other data sets, ph1 . . . phN, where phi=ph0·(fSA0)i+1, for i=1 to N.

9. The method of claim 1, wherein the plurality of other data sets, ph1 . . . phN, comprise a respective plurality of other sub-aperture sets, denoted as SA1 . . . SAN, of the SAR phase history data, PH, and step (c) further comprises multiplying the SAR phase history data, PH, by a respective function representing each one of the plurality of other sub-aperture sets: phi=PH·fSAi, for i=1 to N.

10. The method of claim 8, wherein the plurality of other sub-aperture sets, SA1 . . . SAN, are chosen from one or more of a circularly symmetric Gaussian, a circularly asymmetric Gaussian, a rectangle, an annuli, and a sector.

11. A method according to claim 1, further comprising the steps of: b1) integrating over the first data set, ph0, and each one of the plurality of other data sets, ph1 . . . phN, to obtain a plurality of features, denoted as F0 . . . FN;b2) storing the plurality of features, F0 . . . FN;b3) performing target detection based on the stored plurality of features to obtain a receiver operating characteristic (ROC);b4) determining if the ROC is acceptable; andb5) repeating the method for a new first data set and a new plurality of other data sets if the ROC is unacceptable.

12. The method of claim 11, wherein step (b3) includes inputting the plurality of features to one or more of a statistical, maximum likelihood classifier, and a neural network, packpropagation classifier.

13. The method of claim 12, wherein step (b3) further includes inputting one or more prototype vectors to the classifier(s), each vector corresponding to a respective target or class.

14. The method of claim 11, wherein the first data set, ph0, and the plurality of other data sets, ph1 . . . phN, include a respective plurality of sub-aperture sets, denoted as SA0 . . . SAN, of the SAR phase history data, PH, and the method further includes the step of: choosing the plurality of sub-aperture sets, SA0 . . . SAN, from one or more of a circularly symmetric Gaussian, a circularly asymmetric Gaussian, a rectangle, an annuli, and a sector.

15. A method of generating multi-resolution synthetic aperture radar (SAR) images, the method comprising the steps of: a) obtaining SAR phase history data, denoted as PH;b) computing an inverse Fourier transform of a first data set, denoted as ph0, of the SAR phase history data, PH, to form a highest resolution image, denoted as G0, where G0=FT−1 (ph0);c) successively low-pass filtering the highest resolution image, G0, to form a plurality of correlated lower resolution images, denoted as G1 . . . GN, where Gi+1=LPF(Gi), for i=0 to N−1;d) subtracting successive ones of the highest resolution image, G0, and the lower resolution images, G1 . . . GN, to obtain a plurality of edge enhanced images, denoted as L0 . . . LN−1, where Li=Gi−Gi+1, for i=0 to N−1; ande) obtaining a lowest resolution image, denoted as LN.

16. The method of claim 15, wherein the step of successively low-pass filtering includes successively convolving the highest resolution image, G0, with a weighting function.

17. The method of claim 16, wherein the weighting function is one of a Gaussian weighting function, a triangular weighting function, a rectangular weighting function, a trimodal weighting function, or a broad Gaussian weighting function.

18. A method according to claim 15, further including the step of encoding the plurality of edge enhanced images and the lowest resolution image.

19. The method of claim 18, wherein the encoding step encodes according to one of a Huffman code, range code, arithmetic code, Unary code, Elias code, Fibonacci code, Golomb code, or Rice code.

20. The method of claim 15, wherein the first data set, ph0, comprises the entire SAR phase history data, PH: ph0=PH.

21. The method of claim 15, wherein the first data set, ph0, comprises a first sub-aperture set, denoted as SA0, of the SAR phase history data, PH, and step (b) further includes multiplying the SAR phase history data, PH, by a function representing the first sub-aperture set, denoted as fSA0, where ph0=PH·fSA0.

22. The method of claim 21, wherein the first sub-aperture set, SA0, is chosen from one of a circularly symmetric Gaussian, a circularly asymmetric Gaussian, a rectangle, an annuli, or a sector.

23. A method of iterative target detection in a synthetic aperture radar (SAR) scene, the method comprising the steps of: a) obtaining SAR phase history data, denoted as PH, corresponding to the scene;b) integrating over a plurality of data sets, denoted as ph0 . . . phN, of the phase history data, PH, to obtain a plurality of features, denoted as F0 . . . FN;c) storing the plurality of features, F0 . . . FN;d) performing target detection based on the stored plurality of features to obtain a receiver operating characteristic (ROC) of the scene;e) determining if the ROC is acceptable; andf) repeating steps (b)–(e) for a new plurality of other data sets if the ROC is unacceptable.

24. The method of claim 23, wherein step (b3) includes inputting the plurality of features to one or more of a statistical, maximum likelihood classifier, and a neural network, packpropagation classifier.

25. The method of claim 24, wherein step (b3) further includes inputting one or more prototype vectors to the classifier(s), each vector corresponding to a respective target or class.

26. The method of claim 23, wherein the first data set, ph0, and the plurality of other data sets, ph1 . . . phN, include a respective plurality of sub-aperture sets, denoted as SA0 . . . SAN, of the SAR phase history data, PH, and the method further includes the step of: choosing the plurality of sub-aperture sets, SA0 . . . SAN, from one or more of a circularly symmetric Gaussian, a circularly asymmetric Gaussian, a rectangle, an annuli, and a sector.

27. A method of progressively transmitting multi-resolution synthetic aperture radar (SAR) images, the method comprising the steps of: a) obtaining SAR phase history data corresponding to a scene;b) generating a correlated image pyramid including a plurality of multi-resolution images of the scene;c) generating a Laplacian image pyramid including a plurality of Laplacian images corresponding to a difference of the multi-resolution images in successive levels of the correlated image pyramid;d) encoding each one of the plurality of Laplacian images;e) receiving, from a client, a request for a desired image of the scene;f) transmitting each one of the plurality of Laplacian images, starting from a bottom level of the Laplacian image pyramid, until the desired image of the scene is transmitted.

28. A method according to claim 27, further including the steps of: receiving an indication that the desired image of the scene has been transmitted, wherein step (f) stops transmitting once it receives the indication; andrepeating steps (b)–(f) for a new correlated image pyramid if the indication is not received.

29. A method according to claim 27, further comprising the step of detecting one or more targets in the scene, wherein: step (b) generates the correlated image pyramid comprising a multi-resolution images of the one or more targets;step (e) receives a request for a desired image of a desired target of the one or more targets; andstep (f) transmits until the desired image of the desired target is transmitted.