# Patent

## Title
Method and system for developing and using an image reconstruction algorithm for detecting and imaging moving targets

## Abstract

Image reconstruction approaches that use standard Cartesian-sampled, frequency-domain SAR data, which can be collected by a platform with a single radar antenna, together with its associated metadata as the input and generate an output comprised of a set of images that show the directions of moving targets, estimates of their motions, and focused images of these moving targets.

## Background
The U.S. Government may have a paid-up license in this invention and may have the right in limited circumstances to require the patent owner to license others on reasonable terms as provided for by the terms of a U.S. Government contract.


BACKGROUND OF THE INVENTION
1. Field of the Invention
Generally, the invention involves a process for forming radar images. More specifically, the invention involves a process for detecting moving targets and for forming focused radar images of moving targets.
2. Description of the Related Art
Radar, at its most basic application, is used to measure the range to a target. With knowledge of the speed of propagation of the wave, i.e., electromagnetic wave, that is transmitted toward the target, it is possible to resolve in a first dimension, the distance to the target, based on the received reflected wave or echo. In order to use radar as an imaging tool, it is necessary to collect information about the cross-range of the target, in addition to the first dimension information. This cross-range information is about a second dimension perpendicular to the first dimension.
Synthetic aperture radar (SAR) can be used to collect data in both the first and second dimensions, through a process wherein the reflected waves are measured at different angles with respect to an object-of-interest. This process is referred to in the art as collecting radar measurements over a synthetic (as opposed to a literal) aperture. By taking various measurements of the object-of-interest from varying aspect angles, it is possible to determine approximate distance to the scattering centers within an object-of-interest in the first dimension and location of these scattering centers within the object-of-interest in the second, cross-range dimension. This process of two-dimensional imaging is commonly referred to as reflection tomography.
SAR systems take advantage of the long-range propagation characteristics of radar signals and the complex information processing capability of modern digital electronics to provide high-resolution imagery. SAR imaging is not restricted by time of day or atmospheric conditions due to its operative frequencies. Consequently, SAR imaging supplements other photographic and optical imaging techniques in order to facilitate environmental monitoring, earth-resource mapping, and military operations which may require broad-area imaging at high resolutions. More specifically, SAR technology provides detailed terrain information to geologists for mineral exploration, environmentalists for determination of oil spill boundaries, navigators for sea state and ice hazard mapping, and the military for reconnaissance and targeting information.
Other systems using reflection data, also referred to as projection measurements, are police radars for vehicle speed monitoring and detection, fault inspection systems using acoustic imaging, submarine sonar for imaging underwater objects and the like, seismic imaging system for tunnel detection, oil exploration, geological surveys, etc., and medical diagnostic tools such as sonograms and echocardiograms.
There have been two basic types of processing techniques used in the field of reflection tomography to reconstruct single-bounce (SB) reflection data. First, the frequency-domain projection-slice theorem takes the measured phase history from the reflection data taken at different aspect angles and generates the reconstruction of an image using Fourier transforms. This reconstruction technique is often used for reconstructing SAR image data in order to minimize the computational load that results from necessarily complex processing. A second technique, more prevalent in the medical imaging community, is based on the time-domain back projection techniques. Both of these techniques are discussed in U.S. Pat. No. 5,805,098 to McCorkle which is incorporated herein by reference in its entirety.
An image reconstruction algorithm designed to detect and extract multiple bounce scattering effects on image formation is described in U.S. patent application Ser. No. 10/954,218 entitled PROCESS FOR MAPPING MULTIPLE-BOUNCE GHOSTING ARTIFACTS FROM RADAR IMAGING DATA” filed Oct. 1, 2004 which is a continuation of U.S. patent application Ser. No. 10/631,712 entitled PROCESS FOR MAPPING MULTIPLE-BOUNCE GHOSTING ARTIFACTS FROM RADAR IMAGING DATA” which claims priority to U.S. Pat. No. 6,646,593 similarly titled, which claims priority to U.S. Provisional Patent Application No. 60/345,639, entitled “SPOTLIGHT SAR IMAGE FORMATION WITHOUT MULTIPLE-BOUNCE GHOSTING ARTIFACTS” filed Jan. 8, 2002, all of which are incorporated herein by reference in their entirety.
The reflection data processing techniques described in the related art described herein assume that the impinging wave reflects off of an object of interest that is stationary within the scene of interest before returning back to the receiver. This assumption neglects the situation wherein the wave actually reflects off of a target or object of interest that changes its location within the scene as the sensor is collecting the data used for image reconstruction.
In prior art directed at the task of detecting and imaging moving targets, Perry et al., “SAR Imaging of Moving Targets,” IEEE Aerospace and Electronic Systems, Vol. 35, No. 1, 1999, pp. 188-200, developed a technique for performing SAR imaging of moving targets based upon techniques that are similar to that of conventional moving target indication (MTI) techniques applied to SAR data. In particular, these authors argue that “standard techniques, such as CFAR (constant false alarm rate) detection prescreening may be used to isolate both the static and moving targets from the scene.” They further claim that “this is more difficult for the unfocused movers because they are smeared more than static targets and hence require a larger target-to-clutter ratio for good detection at reasonable false alarm rates.” Thus, these techniques effectively separate the moving targets from the stationary clutter prior to the processing used to actually focus well-formed images of the moving targets. Perry et al. claim “this approach does work, however, for sparse nonurban scenes” as exemplified by their good results for isolated vehicles on a desert road in China Lake, Calif. However, this prior art does not offer such claims of success in more challenging urban environments that typically have a much higher level of stationary clutter. They further claim that the use of “a zero-Doppler clutter filtering technique, such as DPCA (displaced phase center antenna)” can be used to improve the pre-detection of moving targets. However, DPCA requires the platform to have multiple radar antennas.
Other techniques for separating moving targets from SAR data, such as are described in Friedlander et al., “VSAR: A High Resolution Radar System for Detection of Moving Targets,” IEE Proceedings on Radar, Sonar, and Navigation, Vol. 144, No. 4, August 1997, pp. 205-218, are also known. However, these methods require that the platform possess multiple radar antennas.
Other prior art approaches that use conventional SAR data to estimate the locations and velocities of surface targets moving in an arbitrary direction are highly dependent on scene context and, furthermore, require a high level of human interaction to yield viable results. See Carrara et al. “Spotlight Synthetic Aperture Radar-Signal Processing Algorithms,” Artech House, Boston, 1995, pp. 238-241. In particular, the required contextual information is that of a road network. Through use of intensive human analysis and a good SAR simulation tool, it is possible to estimate target location and velocity based on the existence of a cross-range streak in conventional SAR imagery. These approaches use a trial-and-error approach wherein a human analyst performs the following steps: 1) modeling the forward problem, wherein scene context is used to make an initial hypothesis for the location and trajectory of a surface moving target within the scene; 2) using a SAR simulator to analyze the streak signature arising in a conventional SAR image due to the hypothesized moving target; 3) comparing the results of the simulation with the actual measured SAR image; 4) using this comparison to refine the hypothesis of the target location and trajectory; 5) iterating until the human analyst is satisfied with the correlation between the simulated streak signature and the that of the actual measured data; and 6) applying this procedure to the next candidate moving target in the original SAR image, as evidenced by the existence of a cross-range streak in the image. This approach is human-intensive and typically requires auxiliary road network information to yield correct estimates the locations and trajectories of surface moving targets using SAR data. However, such techniques often yield ambiguous results, as there are typically many different roads on which any particular target may be moving and which may give rise to identical signature streaks in the conventional SAR image. Furthermore, certain types of contextual information, e.g. road network information) may not be reliable for many tactical targets, e.g. tanks, because such vehicles are not typically constrained to travel on road networks.
SUMMARY OF THE INVENTION
Summary of the Problem
As discussed above, prior art approaches to detecting and imaging moving targets in conventional SAR imagery suffer various limitations, such as being applicable only to targets moving slowly in the cross-range dimension, requiring multiple antennas, or relying heavily on contextual information and human interaction. These facts motivate the need for the development of a modified image formation process that detects and images moving targets without these limitations.
Summary of the Solution
Described herein are image reconstruction approaches that use standard Cartesian-sampled, frequency-domain SAR data, which can be collected by a platform with a single radar antenna, together with its associated metadata as the input and generate an output comprised of a set of images that show the directions of moving targets, estimates of their motions, and focused images of these moving targets.
In one embodiment of the invention, five basic processing steps implemented together give significant gains over simpler approaches that omit one or more of the steps:
1) Decompose the original Cartesian-sampled, frequency-domain SAR data into some number N of subapertures (possibly overlapping) in the “slow-time” or aperture dimension;
2) Select some number M of candidate target motions vm (m=1, . . . , M) having components in both the down-range radar direction and the cross-range radar direction (and possibly the elevation direction), where these components may be further characterized by one or more parameters (such as velocity, acceleration, etc.);
3) For each value of the N subaperture data sets combined with each value of the candidate target motion vm, form a “motion-focused” subaperture image that is obtained by including additional phase factors into the input Cartesian-sampled, frequency-domain data to obtain a well-focused SAR magnitude image for all surface targets moving with the candidate motion vm;
4) For each value of the candidate motion vm, use the N different motion-focused images as input into an algorithm that jointly computes balanced estimates of the motion-focused subaperture image, which ideally contains only objects moving with motion vm, and the “stationary-defocused” image. This part of the approach utilizes the fact that within the input motion-focused subaperture images, the defocused stationary clutter exhibits a predictable drift with regard to the subaperture index n=1, . . . , N, whereas objects moving with motion vm do not exhibit such drift. This signature difference is exploited to jointly estimate idealized images of the focused objects moving according to motion vm and the defocused stationary clutter; and
5) The resulting subaperture image estimates of the motion-focused targets and the defocused stationary clutter are used to make binary detection decisions regarding whether potential targets moving according to motion vm are present in the image. The value of the assumed candidate motion vm is then used to assign an estimate for the motion of each detected target. The motion-focused subaperture images of the detected targets yield focused images that are similar in appearance to images that would result if these same targets were actually stationary rather than moving within a conventional SAR image.
A further embodiment of the invention is a process for detecting a moving target and reconstructing a scene with a moving target from scattering data, said process comprising: collecting an image of a scene, the image containing scattering data from the scene representing (a) stationary scatterer reflection data and (b) moving scatterer reflection data. The process further includes decomposing said scattering data into a plurality of subsets of scattering data; choosing a plurality of candidate motions; and forming for each pairwise combination of each said subset of scattering data with each said candidate motion a target-focused subimage, wherein said target-focused subimage is focused according to said candidate motion.
For each said candidate motion, the process includes using a plurality of said target-focused subimages focused according to said candidate motion to form joint balanced estimates of (1) a focused subimage containing only scatterers moving with said candidate motion and (2) a defocused image containing only stationary scatterers.
In an embodiment, the data are collected using inverse synthetic aperture radar or related methods wherein the transmitter/receiver locations are stationary during the measurements of the reflected echoes, as with one possible mode of police radars used for vehicle speed detection and monitoring.
In another embodiment, the scattering data are collected based upon the general principles of reflection tomography, as with the monitoring of blood flow using medical ultrasound technology.
In a further embodiment, the estimates of moving entities include two or three components of the entity velocity, as with the estimate of two or three velocity components of blood flow using medical ultrasound technology.
The process further includes using said estimate of the focused subimage containing only scatterers moving with said candidate motion and said estimate of the defocused image containing only stationary scatterers to make a binary decision with regard to the presence or absence of a target moving with said candidate motion.
As one of ordinary skill in the art will recognize, various algorithms may be formulated to implement processing based on the methodology outlined above.



BRIEF DESCRIPTION OF THE DRAWINGS
The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.
In the Figures:
FIG. 1a illustrates true scattering point data from multiple moving and stationary targets according to embodiments of the present invention.
FIG. 1b illustrates a conventional SAR image corresponding to scattering points of FIG. 1a in one embodiment of the present invention.
FIGS. 2a-2b illustrate target-focused images corresponding to scattering points of FIG. 1a according to one embodiment of the present invention.
FIG. 2c illustrates a binary target detection mask corresponding to scattering points of FIG. 1a according to one embodiment of the present invention.
FIGS. 3a-3b illustrate target-focused images corresponding to scattering points of FIG. 1a according to one embodiment of the present invention.
FIG. 4 illustrates a binary target detection mask corresponding to scattering points of FIG. 1a according to one embodiment of the present invention.
FIGS. 5-7 illustrate binary target detection masks corresponding to scattering points of FIG. 1a according to one embodiment of the present invention.



DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS OF THE PRESENT INVENTION
According to an embodiment of the present invention, the detection and imaging of moving targets may be obtained by taking into account candidate target motions through application of an algorithm for moving target estimation via synthetic aperture radar (MESAR). The moving target detection and image reconstruction process begins with complex-valued SAR data in the frequency domain. Such data retains the phase content of the scene reflectivity and is obtained via transmitter and receiver configurations (not shown) that are well-known in the art. The data is processed according to MESAR via a processor (not shown). Processors are well-known in the art and will not be discussed further herein. If the data are provided in the complex-valued image domain, then any existing sidelobe weighting should be removed prior to application of a two-dimensional Fourier transform to obtain Cartesian-resampled frequency domain data (hereinafter “complex Cartesian data”). One of ordinary skill in the art knows how to implement a two-dimensional Fourier transform on complex image data as well as remove sidelobe weighting from that data, hence the details are not included here. In addition, one of ordinary skill in the art also knows how to apply sensor motion compensation and autofocus processes to SAR data, and these processes are assumed to have been applied prior to obtaining the complex Cartesian data.
Let G(ξm,ηn) denote the complex Cartesian data in terms of M discrete resampled values of the down-range spatial frequency ξm and Nf discrete resample values of the cross-range spatial frequency ηn (here the subscript f refers to the original full aperture). The first step in the processing under this embodiment is to use the original full cross-range aperture in order to form a reduced aperture that is smaller than the full aperture. Such a strategy can be beneficial if moving vehicles in the scene exhibit highly non-linear motion over the temporal duration of the full aperture. Let the integer Nr<Nf denote the number of discrete values of the cross-range spatial frequency η for this reduced aperture. Here, the subscript r refers to the reduced aperture.
Given the complex-valued, Cartesian resampled, frequency-domain data Gr(ξm,ηn), which contains Nr samples in the cross-range spatial frequency direction, the approach of this embodiment further selects still smaller subapertures in terms of the cross-range frequency dimension. Select S equally-spaced and overlapping subapertures within each of the reduced subapertures, and let N (where N<Nr) denote the number of samples in the cross-range spatial frequency for each of the S subapertures. This segments the original frequency-domain data G(ξm,ηn), with m=1, . . . , M and n=1, . . . , Nf, into a number of smaller and probably overlapping data sets denoted Gs(ξm,ηn) in terms of the subaperture index s, with m=1, . . . , M and n=1, . . . , N.
Before processing the subaperture data sets, it is necessary to select the candidate target motions, which describe the target motion in both the radar down-range and radar cross-range directions, along with a possible motion component in the elevation direction. For illustration purposes, assume that the target motion in the elevation direction is negligible. One versed in the art could modify the processing below in order to include target motions in the elevation direction, as well.
In the processing under one embodiment of the invention, parameterize the target velocities in terms of time t to obtain the following down-range {tilde over (x)} and cross-range {tilde over (y)} positional components:




  
    
      
        
          
            
              x
              ~
            
            ⁡
            
              (
              t
              )
            
          
          =
          
            
              x
              0
            
            +
            
              
                u
                0
              
              ⁢
              t
            
            +
            
              
                1
                2
              
              ⁢
              
                a
                x
              
              ⁢
              
                t
                2
              
            
          
        
        ,
      
    
    
      
        (
        1
        )
      
    
  
  
    
      
        
          
            y
            ~
          
          ⁡
          
            (
            t
            )
          
        
        =
        
          
            y
            0
          
          +
          
            
              v
              0
            
            ⁢
            t
          
          +
          
            
              1
              2
            
            ⁢
            
              a
              y
            
            ⁢
            
              
                t
                2
              
              .
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




Here the constants x1 and y1 denote the average position of the moving target in the down-range and cross-range directions, respectively, during the coherent processing time of the synthetic aperture. The parameters u0 and v0 determine the vehicle velocity, and the parameters ax and ay determine the acceleration. The time t is assumed to vary between −T0/2 and T0/2. This analysis can easily be extended to include higher order motion effects, such as the t3-order “jerk” motion term. However, such a more complete motion parameterization typically incurs a higher computational burden.

The processing under one embodiment of the invention searches for moving vehicles having some finite number of candidate motion parameter sets characterized by the parameters (uk, vk, ax,k, ay,k), with k=1, . . . , K, with K equal to the total number of candidate velocities. In principle, there is no restriction on the grid of candidate parameter sets. A relatively simple configuration is to define a set of equally-spaced candidate velocities and equally-spaced accelerations on a Cartesian grid in a parameterized velocity-acceleration space. However, it also possible to select candidate motion parameter sets on a non-rectangular grid, for example a grid representing a logarithmic distribution of speeds.
It is necessary to determine the signatures of moving vehicles within SAR images. The goal is to perform this characterization in terms of variables that are convenient for ingesting complex-valued SAR images. To perform this calculation, begin with the following model of the down-converted SAR measurement data:




  
    
      
        
          G
          ⁡
          
            (
            
              f
              ,
              t
            
            )
          
        
        =
        
          
            ∑
            l
          
          ⁢
          
            
              σ
              l
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  
                    -
                    j2π
                  
                  ⁢
                  
                    {
                    
                      
                        R
                        ⁡
                        
                          (
                          
                            
                              θ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            ,
                            
                              φ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                          
                          )
                        
                      
                      -
                      
                        
                          R
                          0
                        
                        ⁡
                        
                          (
                          
                            
                              θ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            ,
                            
                              φ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                          
                          )
                        
                      
                    
                    }
                  
                  ⁢
                  2
                  ⁢
                  
                    f
                    /
                    c
                  
                
                )
              
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




with θ(t) and φ(t) equal to the platform azimuth and elevation angles, respectively, relative to ground plane coordinates. In equation (3), R(θ(t),φ(t)) is equal to the range from the platform to the lth scatterer characterized by complex-valued reflectivity σl, and R0(θ(t),φ(t)) is the range from the platform to the selected ground reference point (GRP) for spotlight SAR image formation. Note that although these equations are presented for spotlight SAR image formation, one versed in the art could modify the equations to accommodate other SAR imaging modes, including stripmap and scan.

The instantaneous difference between R and R0 for the lth scatterer is equal to the projective range from the scatterer to the GRP. This distance can be expressed in the form:

r(θ(t),φ(t))=R(θ(t),φ(t))−R0(θ(t),φ(t))={tilde over (x)}l(t)cos(θ(t))cos(φ(t))+{tilde over (y)}l(t)cos(θ(t))cos(φ(t))  (4)

for wavefronts that are approximately planar at the GRP. It is useful to define a projection of the spatial frequency relative to the projective range via





  
    
      
        ρ
        ≡
        
          
            
              2
              ⁢
              f
            
            c
          
          .
        
      
    
    
      
        (
        5
        )
      
    
  




The down-converted SAR measurements now have the form





  
    
      
        
          G
          ⁡
          
            (
            
              ρ
              ,
              t
            
            )
          
        
        =
        
          
            ∑
            l
          
          ⁢
          
            
              σ
              l
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  
                    -
                    j2π
                  
                  ⁢
                  
                    {
                    
                      
                        
                          
                            
                              x
                              ~
                            
                            l
                          
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                        ⁢
                        
                          cos
                          ⁡
                          
                            (
                            
                              θ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            )
                          
                        
                        ⁢
                        
                          cos
                          ⁡
                          
                            (
                            
                              φ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            )
                          
                        
                      
                      +
                      
                        
                          
                            
                              y
                              ~
                            
                            l
                          
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                        ⁢
                        
                          sin
                          ⁡
                          
                            (
                            
                              θ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            )
                          
                        
                        ⁢
                        
                          cos
                          ⁡
                          
                            (
                            
                              φ
                              ⁡
                              
                                (
                                t
                                )
                              
                            
                            )
                          
                        
                      
                    
                    }
                  
                  ⁢
                  ρ
                
                )
              
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




Equation (6) has a form that facilitates transformation from the spherical spatial frequency domain (ρ,θ,φ) to the Cartesian spatial frequency domain via the transformation:

ξ(t)=ρ cos(θ(t))cos(φ(t)),  (7)

η(t)=ρ cos(θ(t))sin(φ(t)),  (8)

with ξ(t) equal to the down-range spatial frequency, and η(t) equal to the cross-range spatial frequency. Thus, the expression for the idealized SAR phase-history measurements has the form:





  
    
      
        
          
            G
            ⁡
            
              (
              
                f
                ,
                t
              
              )
            
          
          →
          
            
              G
              ~
            
            ⁡
            
              (
              
                ξ
                ,
                η
                ,
                t
              
              )
            
          
        
        =
        
          
            ∑
            l
          
          ⁢
          
            
              σ
              l
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  
                    -
                    j2π
                  
                  ⁢
                  
                    {
                    
                      
                        
                          
                            
                              x
                              ~
                            
                            l
                          
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                        ⁢
                        
                          ξ
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                      
                      +
                      
                        
                          
                            y
                            l
                          
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                        ⁢
                        
                          η
                          ⁡
                          
                            (
                            t
                            )
                          
                        
                      
                    
                    }
                  
                
                )
              
            
          
        
      
    
    
      
        (
        9
        )
      
    
  




The computation of the subaperture SAR images can be performed by evaluating the following expression:




  
    
      
        
          
            b
            s
          
          ⁡
          
            (
            
              x
              ,
              y
            
            )
          
        
        =
        
          
            ∫
            
              
                ξ
                f
              
              -
              
                Δ
                ⁢
                
                  
                
                ⁢
                
                  ξ
                  /
                  2
                
              
            
            
              
                ξ
                f
              
              +
              
                Δξ
                /
                2
              
            
          
          ⁢
          
            
          
          ⁢
          
            
              d
              ξ
            
            ⁢
            
              
                ∫
                
                  
                    η
                    s
                  
                  -
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      η
                      /
                      2
                    
                  
                
                
                  
                    η
                    2
                  
                  +
                  
                    Δη
                    /
                    2
                  
                
              
              ⁢
              
                
              
              ⁢
              
                
                  d
                  η
                
                ⁢
                
                  
                
                ⁢
                
                  
                    G
                    ~
                  
                  ⁡
                  
                    (
                    
                      ξ
                      ,
                      η
                      ,
                      t
                    
                    )
                  
                
                ⁢
                
                  exp
                  ⁡
                  
                    (
                    
                      j2π
                      (
                      
                        
                          x
                          ⁢
                          
                            
                          
                          ⁢
                          ξ
                        
                        +
                        
                          y
                          ⁢
                          
                            
                          
                          ⁢
                          η
                        
                      
                      }
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




with s equal to the subaperture index, and ηs equal to the central cross-range spatial frequency corresponding to subaperture s. A predictive evaluation of equation (10) requires knowledge of the value of the slow-time t for every sample value of the spatial frequency coordinates (ξ,η). In order to perform this calculation, it is useful to parameterize the platform trajectory.

In one embodiment of the invention, a viable parameterization of the platform trajectory is the following in terms of the ground-plane coordinates (X0,Y0,Z0):

X(t)=±V0t sin(θg)cos(λ)−X0,  (11)

Y(t)=±V0t cos(θg)cos(λ),  (12)

Z(t)=V0t sin(λ)+Z0,  (13)

with V0 equal to the platform speed, λ equal to the platform pitch angle relative to level flight, θg equal to the platform ground-plane squint angle, X0 equal to the platform ground range relative to the platform aperture center, and Z0 equal to the platform elevation above the ground plane. The positive sign denotes platform movement in the +y direction, and a minus sign denotes movement in the −y direction.

The aspect angle θ in the spatial domain of the ground-plane coordinates (X0,Y0,Z0) is identical to that of the spatial frequency coordinates (ξ,η,ζ), with ζ equal to the spatial frequency in the elevation direction. Thus the tangent of the aspect angle can be expressed in dual form:




  
    
      
        
          tan
          ⁡
          
            (
            
              θ
              ⁡
              
                (
                t
                )
              
            
            )
          
        
        =
        
          
            
              Y
              ⁡
              
                (
                t
                )
              
            
            
              X
              ⁡
              
                (
                t
                )
              
            
          
          =
          
            
              η
              ξ
            
            .
          
        
      
    
    
      
        (
        14
        )
      
    
  




Equation (14) can be used to solve for the slow-time t in terms of the spatial frequency coordinates (ξ,η) and platform geometrical variables via:




  
    
      
        t
        =
        
          
            ∓
            
              
                X
                0
              
              
                
                  V
                  0
                
                ⁢
                
                  cos
                  ⁡
                  
                    (
                    λ
                    )
                  
                
                ⁢
                
                  cos
                  ⁡
                  
                    (
                    
                      θ
                      g
                    
                    )
                  
                
              
            
          
          ⁢
          
            
              η
              
                {
                
                  ξ
                  -
                  
                    η
                    ⁢
                    
                      
                    
                    ⁢
                    
                      tan
                      ⁡
                      
                        (
                        
                          θ
                          g
                        
                        )
                      
                    
                  
                
                }
              
            
            .
          
        
      
    
    
      
        (
        15
        )
      
    
  




This expression for the slow-time is used in equation (10) to compute subaperture images corresponding to an ideal point target modeled in equations (1) and (2).

The evaluation for the subaperture images bs(x,y) in equation (10) continues in this embodiment by reworking the integral in terms of “centered” spatial frequencies defined by:

{tilde over (ξ)}≡ξ−ξf,  (16)

{tilde over (η)}≡η−ηs.  (17)

The subaperture image equation (10) becomes





  
    
      
        
          
            
              b
              s
            
            ⁡
            
              (
              
                x
                ,
                y
              
              )
            
          
          =
          
            
              σ
              0
            
            ⁢
            
              
                ∫
                
                  -
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    
                      ξ
                      /
                      2
                    
                  
                
                
                  Δξ
                  /
                  2
                
              
              ⁢
              
                
              
              ⁢
              
                
                  d
                  
                    ξ
                    ~
                  
                
                ⁢
                
                  
                    ∫
                    
                      
                        -
                        Δ
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        η
                        /
                        2
                      
                    
                    
                      Δη
                      /
                      2
                    
                  
                  ⁢
                  
                    
                      d
                      
                        η
                        ~
                      
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      exp
                      ⁡
                      
                        (
                        j2πΩ
                        )
                      
                    
                  
                
              
            
          
        
        ,
      
    
    
      
        (
        18
        )
      
    
  




in terms of the phase defined by





  
    
      
        
          Ω
          ≡
          
            
              
                (
                
                  x
                  -
                  
                    x
                    0
                  
                  -
                  
                    
                      u
                      0
                    
                    ⁢
                    t
                  
                  -
                  
                    
                      1
                      2
                    
                    ⁢
                    
                      a
                      x
                    
                    ⁢
                    
                      t
                      2
                    
                  
                
                )
              
              ⁢
              
                (
                
                  
                    ξ
                    f
                  
                  +
                  
                    ξ
                    ~
                  
                
                )
              
            
            +
            
              
                (
                
                  y
                  -
                  
                    y
                    0
                  
                  -
                  
                    
                      v
                      0
                    
                    ⁢
                    t
                  
                  -
                  
                    
                      1
                      2
                    
                    ⁢
                    
                      a
                      y
                    
                    ⁢
                    
                      t
                      2
                    
                  
                
                )
              
              ⁢
              
                (
                
                  
                    η
                    s
                  
                  +
                  
                    η
                    ~
                  
                
                )
              
            
          
        
        ,
      
    
    
      
        (
        19
        )
      
    
  




using equations (1), (2), and (9), with the slow-time t determined by equation (15). For further understanding, it is useful to examine a Taylor series expansion of the terms in equation (18) in terms of three dimensionless parameters that are typically small for narrow-band SAR image formation:





  
    
      
        
          
            ɛ
            x
          
          ≡
          
            
              ξ
              ~
            
            
              ξ
              f
            
          
        
        ,
      
    
    
      
        (
        20
        )
      
    
  
  
    
      
        
          
            ɛ
            y
          
          ≡
          
            
              η
              ~
            
            
              ξ
              f
            
          
        
        ,
      
    
    
      
        (
        21
        )
      
    
  
  
    
      
        
          ɛ
          s
        
        ≡
        
          
            
              η
              s
            
            
              ξ
              f
            
          
          .
        
      
    
    
      
        (
        22
        )
      
    
  




The third parameter is often smaller that the other two for the case of subapertures that are a relatively large fractional size of the full aperture, since there is a relatively small span of possible values for ηs for this case.

Further analysis of the subaperture image equation (18) requires evaluation of the phase in equation (19) through at least second order in the two-dimensional Taylor series expansion in terms of (εx,εy) to assess the signature characteristics of moving scatterers within the SAR imagery. The linear terms will be seen to affect the drift of each mover with regard to subaperture, and the quadratic terms will be seen to affect the focus of each mover. The Taylor series expansion reveals the following for the phase of equation (19):

Ω{x−x0}ξf+{{y−y0}ξf−ku0}εs−k{v0+α0}εs2−kβ0εs3+{{x−x0}ξf+k{v0−α0}εs2}εx+{{y−y0}ξf−ku0−2k{v0+α0}εs−3 kβ0εs2}εy+ku0εsεx2+2k{v0−α0}εsεxεy−k{v0+α0+3β0εs}εy2,  (23)

in terms of the parameters





  
    
      
        
          
            α
            0
          
          ≡
          
            
              
                u
                0
              
              ⁢
              
                tan
                ⁡
                
                  (
                  
                    θ
                    g
                  
                  )
                
              
            
            +
            
              
                k
                ⁢
                
                  
                
                ⁢
                
                  a
                  x
                
              
              
                2
                ⁢
                
                  ξ
                  f
                
              
            
          
        
        ,
      
    
    
      
        (
        24
        )
      
    
  
  
    
      
        
          
            β
            0
          
          ≡
          
            
              
                v
                0
              
              ⁢
              
                tan
                ⁡
                
                  (
                  
                    θ
                    g
                  
                  )
                
              
            
            +
            
              
                k
                ⁢
                
                  
                
                ⁢
                
                  a
                  y
                
              
              
                2
                ⁢
                
                  ξ
                  f
                
              
            
          
        
        ,
      
    
    
      
        (
        25
        )
      
    
  




and the geometry-dependent constant





  
    
      
        k
        ≡
        
          ∓
          
            
              
                
                  X
                  0
                
                ⁢
                
                  ξ
                  f
                
              
              
                
                  V
                  0
                
                ⁢
                
                  cos
                  ⁡
                  
                    (
                    λ
                    )
                  
                
                ⁢
                
                  cos
                  ⁡
                  
                    (
                    
                      θ
                      g
                    
                    )
                  
                
              
            
            .
          
        
      
    
    
      
        (
        26
        )
      
    
  




The maximum intensity of the subaperture SAR images of equation (18) occurs for the position values (x,y) such that the coefficients of εx and εy in equation (23) vanish, giving:




  
    
      
        
          x
          =
          
            
              x
              0
            
            +
            
              
                k
                
                  ξ
                  f
                
              
              ⁢
              
                {
                
                  
                    v
                    0
                  
                  -
                  
                    α
                    0
                  
                
                }
              
              ⁢
              
                ɛ
                s
                2
              
            
          
        
        ,
      
    
    
      
        (
        27
        )
      
    
  
  
    
      
        y
        =
        
          
            y
            0
          
          +
          
            
              k
              
                ξ
                f
              
            
            ⁢
            
              
                {
                
                  
                    u
                    0
                  
                  +
                  
                    2
                    ⁢
                    
                      {
                      
                        
                          v
                          0
                        
                        +
                        
                          α
                          0
                        
                      
                      }
                    
                    ⁢
                    
                      ɛ
                      s
                    
                  
                  +
                  
                    3
                    ⁢
                    
                      β
                      0
                    
                    ⁢
                    
                      ɛ
                      s
                      2
                    
                  
                
                }
              
              .
            
          
        
      
    
    
      
        (
        28
        )
      
    
  




That is, the signature of a moving scatterer characterized by motion parameters (u0,v0,ax,ay) drifts relative to the background of the stationary scatterers with regard to the subaperture center cross-range frequency ηs according to equations (27) and (28). This difference in the drift signature between moving and stationary scatterers may be exploited in the processing of this embodiment of the invention. Of course, equations (27) and (28) can be Taylor expanded through arbitrarily high order in the dimensionless subaperture parameter εs if higher accuracy is required. One versed in the art could perform this expansion to arbitrarily high order.

The processing of this embodiment now searches for targets characterized by one of the candidate motion parameters (uk,vk,ax,k,ay,k), using the S different subapertures. It is first desirable to introduce a focusing term into the frequency-domain data sets of the S subapertures, so as to increase the peak scattering energy of targets that are actually moving with a motion that is determined by the parameters (uk,vk,ax,k,ay,k) within the corresponding subaperture images. An additional benefit of such a quadratic focusing term is that the stationary scatterers become defocused in the resulting subaperture images, so that the peak scattering energy of strong stationary scatterers (e.g. trihedrals) is actually reduced. These two effects contrive to increase the image intensities of targets moving with the candidate motion and to decrease the image intensities of stationary clutter, thereby increasing the effective signal-to-noise ratio of the moving targets of interest
The implementation of the motion-based pre-focusing of the subaperture images is determined by the quadratic phase terms in equation (23) that are proportional to εx2, εxεy, and εy2. The pre-focusing can be accomplished by including additional factors in the subaperture image formation that exactly cancel these quadratic phase terms for the candidate motion parameter set under consideration. Analysis of equation (19) shows that this desired motion-based pre-focusing can be accomplished by including the following factor in the subaperture image formation:




  
    
      
        
          Ω
          ^
        
        ≡
        
          
            exp
            (
            
              
            
            ⁢
            
              j2π
              ⁢
              
                
              
              ⁢
              
                {
                
                  
                    
                      
                        u
                        ^
                      
                      0
                    
                    ⁢
                    t
                    ⁢
                    
                      {
                      
                        
                          ξ
                          f
                        
                        +
                        
                          ξ
                          ~
                        
                      
                      }
                    
                  
                  +
                  
                    
                      
                        v
                        ^
                      
                      0
                    
                    ⁢
                    t
                    ⁢
                    
                      {
                      
                        
                          η
                          s
                        
                        +
                        
                          η
                          ~
                        
                      
                      }
                    
                  
                  +
                  
                    
                      1
                      2
                    
                    ⁢
                    
                      
                        a
                        ^
                      
                      x
                    
                    ⁢
                    
                      t
                      2
                    
                    ⁢
                    
                      {
                      
                        
                          ξ
                          f
                        
                        +
                        
                          ξ
                          ~
                        
                      
                      }
                    
                  
                  +
                  
                    
                      1
                      2
                    
                    ⁢
                    
                      
                        a
                        ^
                      
                      y
                    
                    ⁢
                    
                      t
                      2
                    
                    ⁢
                    
                      {
                      
                        
                          η
                          s
                        
                        +
                        
                          η
                          ~
                        
                      
                      }
                    
                  
                
                }
              
            
            )
          
          .
        
      
    
    
      
        (
        29
        )
      
    
  




with the slow-time t again determined by Equation (15). In addition, it is useful to include the following phase term that is designed to retain the cross-range offset that is due to the down-range velocity:





  
    
      
        
          
            Ω
            ^
          
          offset
        
        ≡
        
          
            exp
            ⁡
            
              (
              
                j2π
                ⁢
                
                  {
                  
                    
                      -
                      
                        
                          k
                          ⁢
                          
                            
                          
                          ⁢
                          
                            
                              u
                              ^
                            
                            0
                          
                        
                        
                          ξ
                          f
                        
                      
                    
                    ⁢
                    
                      {
                      
                        
                          η
                          s
                        
                        +
                        
                          η
                          ~
                        
                      
                      }
                    
                  
                  }
                
              
              )
            
          
          .
        
      
    
    
      
        (
        30
        )
      
    
  




Once û0 is determined accurately and any residual multiple detections of the same moving vehicle are culled, this cross-range offset can be simply corrected by relocating the detection by the cross-range distance





  
    
      
        
          y
          offset
        
        ≡
        
          
            k
            ⁢
            
              
            
            ⁢
            
              
                u
                ^
              
              0
            
          
          
            ξ
            f
          
        
      
    
    
      
        (
        31
        )
      
    
  




to the vehicle's approximate true location.

In summarizing the image refocusing process, the desired focusing in the subaperture images can be obtained by replacing each of the frequency-domain samples of Gs(ξm,ηn) with the following:

Gs,k(ξm,ηn)=Gs(ξm,ηn)exp(j2π{circumflex over (Ω)}),  (32)

in the formation of the focused subaperture images used to search for targets with candidate motion parameters (uk,vk,ax,k,ay,k). Equation (29) has the additional result of offsetting the subaperture drift of scatterers moving with the candidate motion parameters. However, instead each of the defocused stationary scatterers had this same subaperture drift, but in the opposite direction relative to the mover drift in the unfocused subaperture images.

It should be noted that equation (29) can be extended easily to arbitrarily high orders in terms of the parameterized vehicle motion. Consider that the parameterized vehicle motion equations (1) and (2) can be placed in generic form




  
    
      
        
          
            
              x
              ~
            
            ⁡
            
              (
              t
              )
            
          
          =
          
            
              ∑
              
                l
                =
                0
              
              
                L
                x
              
            
            ⁢
            
              
                1
                
                  l
                  !
                
              
              ⁢
              
                μ
                l
              
              ⁢
              
                t
                l
              
            
          
        
        ,
      
    
    
      
        (
        33
        )
      
    
  
  
    
      
        
          
            
              y
              ~
            
            ⁡
            
              (
              t
              )
            
          
          =
          
            
              ∑
              
                m
                =
                0
              
              
                L
                y
              
            
            ⁢
            
              
                1
                
                  m
                  !
                
              
              ⁢
              
                v
                m
              
              ⁢
              
                t
                m
              
            
          
        
        ,
      
    
    
      
        (
        34
        )
      
    
  




with the expansion order integers Lx and Ly permitted to be large but finite. Then the desired focusing factor of Equation (9) is replaced by

{circumflex over (Ω)}0≡exp(j2π{{tilde over (x)}(t){ξf+ξ}+{tilde over (y)}(t){ηs+{tilde over (η)}}}),  (35)

with the slow-time t again determined by Equation (15).

The motion-focused subaperture images have the desired qualities that any target within the scene that is characterized by a given set of candidate motion parameters (uk,vk,ax,k,ay,k) will be well-focused and offset to a location that is invariant with regard to the subaperture index s and lies at the approximate true target location within the scene. An additional benefit is that the scattering energy of the stationary clutter within these target-focused subaperture images is defocused and actually shifts within the imagery in a predictable direction with regard to the subaperture index s, as described below. Thus, the desired moving targets become strongly focused and invariant with regard to subaperture s, whereas the undesired stationary clutter becomes defocused and shifts with regard to subaperture index s. Therefore, the effective signal-to-noise ratio of all moving targets having the candidate motion parameters (uk,vk,ax,k,ay,k) is significantly increased, thereby increasing the probability of correct detection while reducing false alarms due to stationary clutter.
After forming the target-focused subaperture images, the processing of this embodiment utilizes an algorithm that jointly estimates both moving and stationary scattering effects. Specifically, the focused scattering components of targets moving with the candidate motion parameters (uk,vk,ax,k,ay,k) are stationary with regard to the subaperture index s. The defocused scattering component due to stationary clutter can be seen to drift with respect to the subaperture index s according to equations (27) and (28), but with the exact opposite direction. Here, location (x0, y0) equals the mean location of the defocused stationary scattering energy over the extent of the subapertures. The subaperture drift of the defocused stationary clutter is exploited to reduce the effects of false alarms in attempting to detect targets moving with candidate motion parameters (uk,vk,ax,k,ay,k). Consider the difference in the subaperture images between two adjacent subapertures in terms of the index s. Thus, examination of the defocused stationary scattering image intensities for adjacent subaperture images gs+1,k(xm,yn) and gs,k(xm,yn) reveals the following shifts in the down-range and cross-range components, respectively:




  
    
      
        
          
            δ
            x
          
          =
          
            
              -
              
                k
                
                  ξ
                  f
                
              
            
            ⁢
            
              {
              
                
                  v
                  0
                
                -
                
                  α
                  0
                
              
              }
            
            ⁢
            
              
                (
                
                  Δ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ɛ
                    s
                  
                
                )
              
              2
            
          
        
        ,
      
    
    
      
        (
        36
        )
      
    
  
  
    
      
        
          δ
          y
        
        =
        
          
            -
            
              k
              
                ξ
                f
              
            
          
          ⁢
          
            
              {
              
                
                  u
                  0
                
                +
                
                  2
                  ⁢
                  
                    {
                    
                      
                        v
                        0
                      
                      +
                      
                        α
                        0
                      
                    
                    }
                  
                  ⁢
                  
                    (
                    
                      Δɛ
                      s
                    
                    )
                  
                
                +
                
                  3
                  ⁢
                  
                    
                      
                        β
                        0
                      
                      ⁡
                      
                        (
                        
                          Δɛ
                          s
                        
                        )
                      
                    
                    2
                  
                
              
              }
            
            .
          
        
      
    
    
      
        (
        37
        )
      
    
  




Again, the scattering energy of targets moving with motion characterized by (uk,vk,ax,k,ay,k) does not exhibit such drift with regard to the subaperture index s.

To exploit the different subaperture drift signatures of moving targets characterized by (uk,vk,ax,k,ay,k) via the motion index k versus that of the stationary clutter, a measurement model of these effects may be utilized. Let p(xm,yn) denote the ideal image of focused targets moving with motion parameters (uk,vk,ax,k,ay,k), and let q(xm,yn) denote the ideal image of the defocused stationary clutter. Then, one model of measured, target-focused, subaperture images obeys the system of equations:

gs,k(xm,yn)=p(xm,yn)+q(xm+sδx,k,yn+sδy,k)  (38)

with the subaperture index s incremented by unity in progressing from one subaperture to the next adjacent subaperture.

It is possible to extend the signal model of Equation (38) to a more general model that includes deterministic interference due to other movers. Such a more generic model would obey the equations:




  
    
      
        
          
            
              g
              
                s
                ,
                k
              
            
            ⁡
            
              (
              
                
                  x
                  m
                
                ,
                
                  y
                  n
                
              
              )
            
          
          =
          
            
              ∑
              k
            
            ⁢
            
              
                q
                k
              
              ⁡
              
                (
                
                  
                    
                      x
                      m
                    
                    +
                    
                      s
                      ⁢
                      
                        
                      
                      ⁢
                      
                        δ
                        
                          x
                          ,
                          k
                        
                      
                    
                  
                  ,
                  
                    
                      y
                      n
                    
                    +
                    
                      s
                      ⁢
                      
                        
                      
                      ⁢
                      
                        δ
                        
                          y
                          ,
                          k
                        
                      
                    
                  
                
                )
              
            
          
        
        ,
      
    
    
      
        (
        39
        )
      
    
  




in terms of a set of motion fields qk(x,y), with the index k denoting the particular motion field. The field of defocused stationary clutter would be one of these terms in equation (39). One versed in the art could modify the analysis below in order to include these additional effects due to the interference of non-candidate movers when searching for movers having a particular set of candidate motion parameters. However, the majority of the received radar energy is due to stationary scatterers when measuring SAR image data. Thus, equation (38) is sufficient for most applications.

The systems of equations (38) are transformed to the frequency domain so that displacements in the image domain become phase factors in the frequency domain:

Gs,k(ξm,ηn)=Pk(ξm,ηn)+Qk(ξm,ηn)exp(j2π{ξmδs+ηnδy}s).  (40)

The frequency-domain equations (40) are linear in terms of the unknowns (i.e. pk(x,y) and (qk(x,y)) and thus are inverted using standard techniques known to one of ordinary skill in the art to jointly estimate the contributions due to focused moving targets and defocused stationary clutter. This solution is simply based upon the joint pseudo-inverse of the forward problem for the case in which there are a greater number of subaperture images that there are motion models to be estimated (i.e. the case in which the number S of subaperture images in greater than 2).

The resulting frequency domain estimates of the images corresponding to targets moving with candidate motion parameters (uk,vk,ax,k,ay,k) and of the stationary clutter are jointly estimated to be the following:




  
    
      
        
          
            
              
                P
                ^
              
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          =
          
            
              
                
                  
                    A
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
                ⁢
                S
              
              -
              
                
                  
                    B
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
                ⁢
                
                  C
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
              
            
            
              
                S
                2
              
              -
              
                
                  C
                  2
                
                ⁡
                
                  (
                  
                    
                      ξ
                      m
                    
                    ,
                    
                      η
                      n
                    
                  
                  )
                
              
            
          
        
        ,
      
    
    
      
        (
        41
        )
      
    
  
  
    
      
        
          
            
              
                Q
                ^
              
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          =
          
            
              
                
                  -
                  
                    
                      A
                      k
                    
                    ⁡
                    
                      (
                      
                        
                          ξ
                          m
                        
                        ,
                        
                          η
                          n
                        
                      
                      )
                    
                  
                
                ⁢
                
                  C
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
              
              +
              
                
                  
                    B
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
                ⁢
                S
              
            
            
              
                S
                2
              
              -
              
                
                  C
                  2
                
                ⁡
                
                  (
                  
                    
                      ξ
                      m
                    
                    ,
                    
                      η
                      n
                    
                  
                  )
                
              
            
          
        
        ,
      
    
    
      
        (
        42
        )
      
    
  




in terms of the functions





  
    
      
        
          
            
              A
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          ≡
          
            
              ∑
              s
            
            ⁢
            
              
                G
                
                  s
                  ,
                  k
                
              
              ⁡
              
                (
                
                  
                    ξ
                    m
                  
                  ,
                  
                    η
                    n
                  
                
                )
              
            
          
        
        ,
      
    
    
      
        (
        43
        )
      
    
  
  
    
      
        
          
            
              B
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          ≡
          
            
              ∑
              s
            
            ⁢
            
              
                
                  G
                  
                    s
                    ,
                    k
                  
                
                ⁡
                
                  (
                  
                    
                      ξ
                      m
                    
                    ,
                    
                      η
                      n
                    
                  
                  )
                
              
              ⁢
              
                exp
                ⁡
                
                  (
                  
                    
                      -
                      j2π
                    
                    ⁢
                    
                      {
                      
                        
                          
                            ξ
                            m
                          
                          ⁢
                          
                            δ
                            x
                          
                        
                        +
                        
                          
                            η
                            n
                          
                          ⁢
                          
                            δ
                            y
                          
                        
                      
                      }
                    
                    ⁢
                    s
                  
                  )
                
              
            
          
        
        ,
      
    
    
      
        (
        44
        )
      
    
  
  
    
      
        
          
            C
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          ≡
          
            
              ∑
              s
            
            ⁢
            
              exp
              ⁡
              
                (
                
                  j2π
                  ⁢
                  
                    {
                    
                      
                        
                          ξ
                          m
                        
                        ⁢
                        
                          δ
                          x
                        
                      
                      +
                      
                        
                          η
                          n
                        
                        ⁢
                        
                          δ
                          y
                        
                      
                    
                    }
                  
                  ⁢
                  s
                
                )
              
            
          
        
        ,
      
    
    
      
        (
        45
        )
      
    
  




with the summations in Eqs. (43)-(45) applied over the subaperture index s and the constant S denoting the total number of subaperture data sets used in these summations.

Finally, the two-dimensional inverse Fourier transforms of equations (41) and (42) are used to obtain the joint estimate of the focused moving targets pk(xm,yn) and the defocused stationary clutter qk(xm,yn):

{circumflex over (p)}k(xm,yn){circumflex over (P)}k(ξm,ηn),  (46)

{circumflex over (q)}k(xm,yn){circumflex over (Q)}k(ξm,ηn).  (47)

The solution for pk(xm,yn) and qk(xm,yn) as determined by equations (41)-(47) is not well-defined for spatial frequencies such that ξδs+ηδy equals an integer, as shown by the denominators of equations (41) and (42) via equation (45). Well-defined solutions exist for spatial frequencies such that ξδs+ηδy is not equal to an integer; however, this problem becomes ill-conditioned for some non-zero span of frequencies in the neighborhood of the locus where ξδs+ηδy equals an integer for the realistic case of noisy estimates of the various subaperture images gs,k(xm,yn) resulting from noisy radar echo measurements. Thus, this solution tends to be sensitive to additive noise in the measured data.
One practical means for ameliorating this ill-conditioning problem is to apply an appropriate weighting function Wk(ξm,ηn) to Pk(ξm,ηn) and Qk(ξm,ηn) in equations (41) and (42) prior to applying the inverse Fourier transform of equations (46) and (47) to obtain the desired functions pk(xm,yn) and qk(xm,yn). Empirical analysis has revealed that the application of a two-dimensional “zero-forcing” window defined by




  
    
      
        
          
            W
            k
          
          ⁡
          
            (
            
              
                ξ
                m
              
              ,
              
                η
                n
              
            
            )
          
        
        ≡
        
          {
          
            
              
                
                  
                    1
                    ,
                    
                      
                        
                          if
                          ⁢
                          
                            
                          
                          ⁢
                          
                            S
                            2
                          
                        
                        -
                        
                          
                            C
                            2
                          
                          ⁡
                          
                            (
                            
                              
                                ξ
                                m
                              
                              ,
                              
                                η
                                n
                              
                            
                            )
                          
                        
                      
                      ≥
                      
                        W
                        min
                      
                    
                  
                
              
              
                
                  
                    0
                    ,
                    
                      
                        
                          if
                          ⁢
                          
                            
                          
                          ⁢
                          
                            S
                            2
                          
                        
                        -
                        
                          
                            C
                            2
                          
                          ⁡
                          
                            (
                            
                              
                                ξ
                                m
                              
                              ,
                              
                                η
                                n
                              
                            
                            )
                          
                        
                      
                      <
                      
                        W
                        min
                      
                    
                  
                
              
            
            ,
          
        
      
    
    
      
        (
        48
        )
      
    
  




with Wmin equal to some appropriate value. An analysis based on the signal-to-noise ratio indicates that the value of Wmin should be set so that the denominator S2−C2(ξm,ηn) of equations (41) and (42) does not vary over a greater span than does the signal-to-noise ratio. Thus, the value of Wmin can be set according to the following:





  
    
      
        
          
            W
            min
          
          ≡
          
            
              max
              ⁢
              
                {
                
                  
                    W
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        ξ
                        m
                      
                      ,
                      
                        η
                        n
                      
                    
                    )
                  
                
                }
              
            
            ψ
          
        
        ,
      
    
    
      
        (
        49
        )
      
    
  




with ψ equal to the estimate of the peak signal-to-noise ratio in absolute (i.e. non-decibel) units. Further analysis has revealed that additional conditioning can be obtained by the application of a two-dimensional Hanning window

Wk,Hann(ξm,ηn)≡{4/3}{1−cos(2π{ξmδx+ηnδy})}  (50)

over the entire span of both ξm and ηn. In summary, one possible conditioning window is

Wk(ξm,ηn)≡Wk,zero(ξm,ηn)Wk,Hann(ξm,ηn)  (51)

in terms of Equations (49) and (50). One versed in the art can examine the application of any number of known techniques for conditioning the system of equations given by Equations (41)-(47).

The ill-conditioning of the original system of equations can be mitigated by using the following estimates for the Fourier transforms of the desired functions pk(xm,yn) and qk(xm,yn), respectively:




  
    
      
        
          
            
              
                P
                ^
              
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          =
          
            
              
                W
                k
              
              ⁡
              
                (
                
                  
                    ξ
                    m
                  
                  ,
                  
                    η
                    n
                  
                
                )
              
            
            ⁢
            
              {
              
                
                  
                    
                      
                        A
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            ξ
                            m
                          
                          ,
                          
                            η
                            n
                          
                        
                        )
                      
                    
                    ⁢
                    S
                  
                  -
                  
                    
                      
                        B
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            ξ
                            m
                          
                          ,
                          
                            η
                            n
                          
                        
                        )
                      
                    
                    ⁢
                    
                      C
                      ⁡
                      
                        (
                        
                          
                            ξ
                            m
                          
                          ,
                          
                            η
                            n
                          
                        
                        )
                      
                    
                  
                
                
                  
                    S
                    2
                  
                  -
                  
                    
                      C
                      2
                    
                    ⁡
                    
                      (
                      
                        
                          ξ
                          m
                        
                        ,
                        
                          η
                          n
                        
                      
                      )
                    
                  
                
              
              }
            
          
        
        ,
      
    
    
      
        (
        52
        )
      
    
  
  
    
      
        
          
            
              Q
              ^
            
            k
          
          ⁡
          
            (
            
              
                ξ
                m
              
              ,
              
                η
                n
              
            
            )
          
        
        =
        
          
            
              W
              k
            
            ⁡
            
              (
              
                
                  ξ
                  m
                
                ,
                
                  η
                  n
                
              
              )
            
          
          ⁢
          
            
              {
              
                
                  
                    
                      -
                      
                        
                          A
                          k
                        
                        ⁡
                        
                          (
                          
                            
                              ξ
                              m
                            
                            ,
                            
                              η
                              n
                            
                          
                          )
                        
                      
                    
                    ⁢
                    
                      C
                      ⁡
                      
                        (
                        
                          
                            ξ
                            m
                          
                          ,
                          
                            η
                            n
                          
                        
                        )
                      
                    
                  
                  +
                  
                    
                      
                        B
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            ξ
                            m
                          
                          ,
                          
                            η
                            n
                          
                        
                        )
                      
                    
                    ⁢
                    S
                  
                
                
                  
                    S
                    2
                  
                  -
                  
                    
                      C
                      2
                    
                    ⁡
                    
                      (
                      
                        
                          ξ
                          m
                        
                        ,
                        
                          η
                          n
                        
                      
                      )
                    
                  
                
              
              }
            
            .
          
        
      
    
    
      
        (
        53
        )
      
    
  




The primary cost in using a weighting function for system conditioning is that the resulting estimates for the focused moving targets pk(xm,yn) and the stationary clutter qk(xm,yn) exhibit some loss in resolution. However, this slight loss of resolution is of little importance to the problem of detecting and imaging moving targets.

In attempting to separate the image of focused targets moving with motion index k from that of the defocused stationary clutter, it is useful to compute the following discrimination image:




  
    
      
        
          
            Λ
            k
          
          ⁡
          
            (
            
              
                x
                m
              
              ,
              
                y
                n
              
            
            )
          
        
        ≡
        
          
            
              
                
                
                  
                    
                      p
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        x
                        m
                      
                      ,
                      
                        y
                        n
                      
                    
                    )
                  
                
                
              
              2
            
            
              
                
                
                  
                    
                      q
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        x
                        m
                      
                      ,
                      
                        y
                        n
                      
                    
                    )
                  
                
                
              
              2
            
          
          .
        
      
    
    
      
        (
        54
        )
      
    
  




This function differentiates the effects of focused targets moving with motion index k from that of the defocused stationary clutter. In particular, the image Λk(xm,yn) is expected to be greater than some constant Λmin>1 for spatial pixels for which the scattering intensities of the k-motion parameter targets are stronger than that due to the defocused stationary clutter, and this image is expected to be less than Λmin for pixels with stronger stationary scattering clutter.

To make a binary declaration (i.e. target present or target not present) of a moving target with motion index k, one approach under this embodiment of the invention is to require Λk(xm,yn) to be greater than Λmin and that pk(xm,yn) be greater than some minimum energy threshold. This processing yields a binary target detection image for motion index k:

Γk(xm,yn)≡{Λk(xm,yn)>Λmin}&{|{circumflex over (p)}k(xm,yn)|>pmin}.  (55)

This binary detection image can be used as an overlay on the original input image.

To obtain images of the actual moving targets, it is possible to use the pk(xm,yn) computed earlier in the processing to obtain an image overlay with slightly degraded resolution. As an alternative, it is possible use the estimates of the candidate motion parameters corresponding to the highest value of the MESAR discrimination image in equation (54) into the quadratic focus factors of equation (29) applied to the original full aperture data. Such processing can be used in order to obtain an image overlay of the moving targets with no degradation in resolution.
It should also be clarified that Equation (41) can be expressed in the generic form of Wiener filtering:

G=ΣAQ+N.  (56)

It is well-known to those versed in the art that such Wiener filtering techniques can offer improved results if the noise spectrum is non-uniform. However, in applying Wiener filtering for the current application, proper mitigation of the ill-conditioning due to the deterministic clutter is still required (e.g., Equation (51)). In addition, the increased computation time of applying Wiener filtering must be traded against that required in order to investigate different motion hypotheses.

The moving target detection and image reconstruction processes developed above are applied to simulated phase history data that contains contributions from both stationary scattering points and moving scattering points with the following results. FIG. 1a is a “truth” plot 100 showing an example containing eighteen stationary point scatterers, as denoted by the red circles 101. The location of three moving point scatterers are denoted by the curved lines, which begin at the green x's at the beginning of the simulate SAR measurements and end at the blue asterisks at the end of the measurements. Thus, the lengths of these lines indicates the true target speeds. The simulation represented in FIG. 1a includes point targets moving at various speeds and in various directions. The image powers of the stationary points were all chosen to be approximately 1 dB higher than that of the moving point targets. This presents a challenging scenario in which stationary clutter can interfere with the detection of the moving targets. The specific values of the truth target velocity and acceleration parameters (uk,vk,ax,k,ay,k) are given by (2, −10, 0.1, −0.2), (4, −2, −0.2, −0.1), and (−2, −6, 0.1, 0.1). In the analysis here and below, all velocities are given in meters/second and all accelerations are in meters/second/second. The components of the velocity and acceleration in the elevation direction (i.e., the z direction) are zero.
The true values of the moving and stationary scattering centers were used to simulate the phase history measurements for a signal-to-noise ratio of 50 dB. The depression angle of the simulated collection was set equal to zero, and the squint angle was also set equal to zero. The simulated waveforms covered a frequency span of 9.75 GHz through 10.25 GHz, and the simulated aperture angles covered a span of −3 degrees through +3 degrees. To complete the preprocessing, these phase history measurements were resampled onto a Cartesian frequency-domain grid to yield frequency-domain data G(ξm,ηn), with m=1, . . . , M, and n=1, . . . , Nf, which is the input data into the overall processing. These simulated data yield the conventional SAR image 102 shown in FIG. 1b, which shows that the stationary targets yield well-focused points and the moving targets yield various smeared energy contributions.
The processing is initiated based upon the full aperture frequency-domain data G(ξm,ηn) m=1, . . . , M, and n=1, . . . , Nf. The original full aperture is decomposed into a total of S=16 equally-spaced and overlapping subapertures, each having 50% of the sample points in the cross-range spatial frequency (i.e., N=Nf/2). A number of candidate sets of target velocity and acceleration parameters (uk,vk,ax,k,ay,k) is selected for analysis. First, we examine the results obtained when the candidate motion parameter set is not matched to the truth. FIG. 2a and FIG. 2b show examples of the first and last (#16) target-focused subaperture images corresponding to the motion parameters (uk,vk,ax,k,ay,k)=(−2,6,−0.1,−0.2). FIGS. 2a and 2b show the expected relative shift in the defocused stationary clutter. FIG. 2c shows that the resulting binary mask Γk(xm,yn) of equation (55) correctly yields no target detections of moving targets characterized by this set of motion parameters.
This same processing is applied to search for moving targets characterized by various values of the candidate motion parameters. FIGS. 3a, 3b, and 4 examine the particular case of (uk,vk,ax,k,ay,k)=(2,−20,0.1,−0.2), which corresponds to one of the moving targets in the original simulated data. FIGS. 3a and 3b reveal the shift of the defocused stationary scattering energy in the first and last (#16) target-focused subaperture images. FIG. 4 shows that the resulting binary mask Γk(xm,yn) for this parameter set correctly detects the true moving target. FIGS. 5-6 show the resulting binary masks for the other two correctly detected targets. FIG. 7 shows the final combined moving target detection mask images resulting from the binary “and” operation applied to all of the binary masks Γk(xm,yn) for the different motion parameter sets (uk,vk,ax,k,ay,k). Thus, this simulation consisting of both moving targets and stationary clutter detected all three moving targets with no false alarms, for one particular setting of the energy threshold pmin in equation (55).
The specific examples described above are not intended to be limiting. Alteration of the processes described above is contemplated by one skilled in the art for application of the MESAR processing to various types of reflection tomography data, including radar, sonar, and the like, for embodiments in which the transmitter/receiver locations move during the measurements of the reflected echoes, as with SAR, or the transmitter/receiver locations do not move, as with Inverse SAR (ISAR).

## Claims

1. A process for detecting a moving target and reconstructing a scene with a moving target from scattering data, said process comprising:
collecting an image of a scene, the image containing scattering data from the scene representing
(a) stationary scatterer reflection data and
(b) moving scatterer reflection data;

decomposing said scattering data into a plurality of subsets of scattering data;
choosing a plurality of candidate motions;
forming for each pairwise combination of each said subset of scattering data with each said candidate motion a target-focused subimage, wherein said target-focused subimage is focused according to said candidate motion;
for each said candidate motion, using a plurality of said target-focused subimages focused according to said candidate motion to form joint balanced estimates of (1) a focused subimage containing only scatterers moving with said candidate motion and (2) a defocused image containing only stationary scatterers;
for each said candidate motion using said estimate of the focused subimage containing only scatterers moving with said candidate motion and said estimate of the defocused image containing only stationary scatterers to make a binary decision with regard to the presence or absence of a target moving with said candidate motion.

2. The process according to claim 1, wherein the scattering data is two-dimensional, Cartesian-sampled, complex-valued, spatial frequency domain data.

3. The process according to claim 2, wherein one of the two dimensions is a cross-range dimension, and decomposing said scattering data into a plurality of subsets of scattering data is in the cross-range dimension.

4. The process according to claim 1, wherein each candidate motion comprises a plurality of components.

5. The process according to claim 1, wherein the image is collected using synthetic aperture radar.

6. The process according to claim 1, wherein the scattering data are collected using a co-located transmitter and receiver.

7. The process according to claim 1, wherein the scattering data are collected using inverse synthetic aperture radar.

8. The process according to claim 1, wherein the scattering data are collected using a method wherein transmitter or receiver locations are stationary during a measurement of a reflected echo.

9. The process according to claim 1, wherein the scattering data are collected for vehicle speed detection or monitoring.

10. The process according to claim 1, wherein the scattering data are collected based upon the general principles of reflection tomography.

11. The process according to claim 1, wherein the scattering data are collected for monitoring of blood flow using medical ultrasound technology.

12. The process according to claim 1, wherein an estimate of a moving entity includes two or more components of a velocity of the entity.

13. The process according to claim 12, wherein the estimate is of blood flow using medical ultrasound technology.

14. A system for detecting a moving target and reconstructing a scene with a moving target from scattering data, said system comprising:
a radar for collecting an image of a scene, the image containing scattering data from the scene representing
(a) stationary scatterer reflection data and
(b) moving scatterer reflection data;

wherein said scattering data is decomposed into a plurality of subsets of scattering data and a plurality of candidate motions is chosen;
for each pairwise combination of each said subset of scattering data with each said candidate motion forming a target-focused subimage, wherein said target-focused subimage is focused according to said candidate motion;
for each said candidate motion, using a plurality of said target-focused subimages focused according to said candidate motion to form joint balanced estimates of (1) a focused subimage containing only scatterers moving with said candidate motion and (2) a defocused image containing only stationary scatterers;
for each said candidate motion using said estimate of the focused subimage containing only scatterers moving with said candidate motion and said estimate of the defocused image containing only stationary scatterers to make a binary decision with regard to the presence or absence of a target moving with said candidate motion.

15. The system according to claim 14 wherein the scattering data is two-dimensional, Cartesian-sampled, complex-valued, spatial frequency domain data.

16. The system according to claim 15, wherein one of the two dimensions is a cross-range dimension, and decomposing said scattering data into a plurality of subsets of scattering data is in the cross-range dimension.

17. The system according to claim 14, wherein each candidate motion comprises a plurality of components.

18. The system according to claim 14, wherein the radar is a synthetic aperture radar.

19. The system according to claim 14, wherein the scattering data is collected using a co-located transmitter and receiver.

20. The system according to claim 14, wherein the radar is an inverse synthetic aperture radar.