Abstract:
In the present invention, the histogram model used in H-PMHT is extended to treat the problem of tracking using hyper-spectral data. Completely general spectral density functions are handled via the use of non-parametric methods. The present invention is not restricted to derivations based on knowledge of the spectral character of the source being tracked. The source spectrum can be estimated in a non-parametric fashion based on an initial track, and this allows the invention to adapt to the source spectrum in situ. The resulting method has improved crossing track performance on sources that have some degree of spectral distinction and will perform no worse than regular H-PMHT on sources that have identical spectral densities.

Description:
STATEMENT OF GOVERNMENT INTEREST 
   The invention described herein may be manufactured and used by or for the Government of the United States of America for governmental purposes without the payment of any royalties thereon or therefore. 
   CROSS REFERENCE TO OTHER PATENT APPLICATIONS 
   Not applicable. 
   BACKGROUND OF THE INVENTION 
   (1) Field of the Invention 
   The present invention relates to remote sensing and remote imaging, and more specifically to a method for processing hyper-spectral remote sensor data for the purpose of displaying the spatial tracks of energy sources in multi-spectral images corresponding to the sensor data. 
   (2) Description of the Prior Art 
   Remote sensing of the energy signals of a moving vehicle or energy source for the purposes of tracking the vehicle or energy source has often been accomplished by measuring the intensity of the energy signals with sensors specifically designed to detect energy intensity. In remote sensing applications the sensor is often a planar array of sensing cells, each cell responding to the energy incident on its corresponding section of the array surface. In other applications, such as acoustic sensing, the received energy on sensor elements must be interpreted through a beam forming function to yield energy intensity in a set of spatially directed cells (more commonly called beams). Such a method is designed to track energy peaks as they move over time on the given set of sensor cells. The total broadband energy is plotted and visually displayed. Targets appear as peaks of energy in the display, and are tracked. One method of target tracking based on sensor level data is the Histogram Probabilistic Multi-Hypothesis Tracking (H-PMHT) algorithm. It is an application of the Expectation-Maximization (EM) method of target tracking. It uses a synthetic (multi-dimensional) histogram interpretation of the received power levels in all of the sensor cells. The data for the H-PMHT algorithm usually consists of broadband intensities on a set of spatial sensor cells. The H-PMHT algorithm has its limitations. For example, in situations where more than one target is being tracked and the targets cross paths, the intensity of the energy signals of the targets merge, making it impossible to distinguish the energy between the two targets. In such a situation, the targets must be reacquired by the sensors after they have crossed resulting in a gap and delay in tracking information. 
   U.S. patent application Ser. No. 10/214,551 to Struzinski teaches a method and system for predicting and detecting the crossing of two target tracks in a bearing versus time coordinate frame. The method/system uses a series of periodic bearing measurements of the two target tracks to determine a bearing rate and a projected intercept with a bearing axis of the bearing versus time coordinate frame. A crossing time t c  for the two target tracks is determined using the tracks&#39; bearing rates and projected intercepts. A prediction that the two target tracks will cross results if a first inequality is satisfied while a detection that the two target tracks have crossed results if a second inequality is satisfied. This method does not, however, address the problem of distinguishing between and identifying the targets before, during and after they have crossed. 
   There is currently no reliable method by which targets can be consistently tracked and distinguished as they cross paths. What is needed is a method for tracking targets that does not rely solely on broadband energy signal intensity, but also utilizes the spectral aspects of the energy signal, combining both intensity and spectral data so that crossing targets can be tracked provided they have some degree of spectral distinction. 
   SUMMARY OF THE INVENTION 
   It is a general purpose and object of the present invention to provide a method for tracking both the spatial sensor data and hyper-spectral sensor data associated with a target. 
   It is a further purpose to estimate a frequency spectrum for a target contribution that consists of only the target&#39;s energy contribution as opposed to the target energy and the noise energy together. 
   These objects are accomplished with the present invention by taking the histogram model used in H-PMHT and extending it to treat the problem of tracking using hyper-spectral data. In the present invention each measurement scan is now a multi-dimensional array wherein each spatial cell has an associated vector of amplitudes in several (possibly disjoint) spectral cells. The intensity data in the multi-dimensional array is interpreted as the spatial-spectral histogram of a synthetic shot process. A statistical model of the random variation of individual cell intensities from scan to scan is required. The procedure adopted in H-PMHT is to quantize the data vector into a “pseudo-histogram,” and then use a multinomial distribution to model the cell counts where the PMHT target mixture model parameterizes the multinomial distribution. The target mixture model determines the cell probabilities that correspond to expected cell counts. 
   The present invention modifies H-PMHT by using a non-parametric spectral characterization of the energy intensity of the target that is assumed known. The use of such a spectral template enhances low signal to noise ratio (SNR) tracking and allows discrimination of spectrally distinct sources as they cross in the spatial domain. The track solutions from previous batches are used to estimate the (non-parametric) spectral characterization that is used to initiate the generation of the updated solutions as new data is received and processed. 
   The present invention provides a mechanism for separating the observed hyperspectral energy into the hyperspectral energy for each source (including noise) using the known spectral characteristics for each source. Completely general spectral density functions are handled via the use of non-parametric methods. In the alternative, the source spectrum is estimated in a non-parametric fashion based on an initial track, allowing the algorithm to adapt to the source spectrum in situ. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     A more complete understanding of the invention and many of the attendant advantages thereto will be readily appreciated as the same becomes better understood by reference to the following detailed description when considered in conjunction with the accompanying drawings depicting an underwater application as the preferred embodiment wherein: 
       FIG. 1  shows an underwater vehicle towing sensors; 
       FIG. 2  shows sensors detecting energy from different sources; 
       FIG. 3  shows the data cube created after raw sensor data is processed; and 
       FIG. 4  shows a flow chart of the method. 
   

   DESCRIPTION OF THE PREFERRED EMBODIMENT 
   Referring now to  FIG. 1  there is shown an underwater vehicle  10  towing an array of sensors  15  arranged on a cable  20 . The sensors  15  are of a type known by those skilled in the art of signal processing such as hydrophones. The sensors  15  are capable of detecting energy signals and their intensities from different directions as illustrated in  FIG. 2 , which shows two vessels labeled k 1  and k 2  and the energy signals  17  emanating from the vessels. The sensor data from each sensor  15  is transmitted along cable  20  to data processors (not shown) within underwater vehicle  10 . The data processors take the raw sensor data and create a data cube  25  as illustrated in  FIG. 3 . Each such data cube  25  is a collection of smaller cubes referred to as sensor cells  30  which correspond to the processed sensor data generated by the sensors  15 . Each sensor cell  30  contains spatial measurements along the x-axis  31 , spectral measurements along the y-axis  32  and time measurements along the t-axis  33 . The side of each sensor cell  30  contained in the (x,t) plane corresponding to spatial measurement is referred to as a spatial cell  36 . The side of each sensor cell  30  contained in the (y,t) plane corresponding to spectral measurement is referred to as a spectral cell  38 . The processing and arrangement of the raw sensor data into a data cube  25  composed of multiple sensor cells  30  that are further composed of spatial cells  36  and spectral cells  38  is known by those skilled in the art of signal processing and is achieved by what is often termed beamforming followed by spectral analysis of the beam intensity data using standard discrete Fourier transform (DFT) techniques known in the art. A single layer of the data cube  25  is referred to as a scan of the sensor space  35  as illustrated in  FIG. 3 . 
   In the preferred embodiment of the present invention, let C={C 1 , . . . , C s }, S≧1, denote the collection of all possible sensor cells  30 . It is assumed that C i ∩C j =φ for all i and j and that C 1 ∪ . . . ∪C s =R dim(c) , where dim(C) denotes the dimension of the sensor space  35 . Furthermore, the sensor cells  30  C={C 1 , . . . , C s } are the Cartesian products of U disjoint spatial cells  36  {D 1  . . . , D U } and V disjoint spectral cells  38  {ε 1 , . . . , ε V }. This particular choice of spatial  36  and spectral  38  cells is application dependent, but they are intrinsically fixed. The total number of sensor cells  30  in a scan of the sensor space  35  is S=UV, and every cell C 1  can be written in the form
 
 C   l   =D   i ×ε j  
 
for some (unique) choice of the cells D i  and ε j . Let D=D 1 ∪ . . . ∪D u  and ε=ε 1 ∪ . . . ∪ε V . The sensor cells  30  from which measurements are available may vary from scan to scan. The sensor cells  30  displayed at time t, for the scan of the sensor space  35 , are the Cartesian product of the spatial cells  36  {D 1 (t), . . . , D U(t) (t)} and the spectral cells  38  {E 1 (t), . . . , E V(t) (t)}, so that the (i,j) th  sensor cell  30  is C ij (t)=D i (t)×E j (t), where 1≦U(t)≦U and 1≦V(t)≦V. The remaining sensor cells  30  in the scan of sensor space  35  are said to be truncated, and no measurements are collected for these cells at time t.
 
   Let the scan of sensor space  35  at time t be denoted by
 
 Z   t   ={z   t,l,1   , . . . , z   t,U(t),V(t) },
 
where z tij≧0  is the output of the sensor space  35  at time t in cell C ij (t), i=1, 2, . . . U(t), j=1, 2, . . . V(t). Let  h   2 &gt;0 be a specified quantization level, and let n tij  denote the quantized value corresponding to the intensity z tij  in cell C ij (t), where
 
                     n   tij     =     ⌊       z   tij       h   2       ⌋       ,           (   1   )               
and |_x_| denotes the greatest integer less than or equal to x. The use of the quantized values {n tij } instead of the measurements {z tij } is an intermediate step in the development. After deriving the auxiliary function of the H-PMHT algorithm using the synthetic counts {n tij }, the measurements {z tij } are recovered in the limit as  h   2 →0.
 
   The “rectangular” spatial-spectral sensor cell structure enables simplifications to the basic equations of H-PMHT. These equations are restated here with the updated notation corresponding to this new cell structure. The cell probability, P ij , for the (i,j) th  cell  30  takes the form 
                       P   ij     ⁡     (     X   t     )       =       ∫     C   ij     (   t   )         ⁢       f   ⁡     (     u   ,     v   ❘     X   t         )       ⁢           ⁢     ⅆ   u           ,     ⅆ   v     ,           (   2   )               
where the sample Probability Density Function (PDF) f(u, v|X t ) is defined over all (u, v)εR dimD ×R dimε =R dimC , by the mixture density
 
                   f   ⁡     (     u   ,     v   ❘     X   t         )       =       ∑     k   =   0     M     ⁢       π   tk     ⁢       G   k     ⁡     (     u   ,     v   ❘     x   tk         )                   (   3   )               
and where π tk  is the component mixing proportion, X t ={x t0 , . . . , x tM } are the component spatial state parameters and G k (u, v|x tk ) is the component PDF corresponding to target k if k≧1 and to noise if k=0. The expected sensor space measurement  z   tij  takes the form
 
                   z   tij     =     {           z   tij           {             1   ≤   i   ≤     U   ⁡     (   t   )         ,                 1   ≤   j   ≤     V   ⁡     (   t   )         ,                              z   t          ⁢         P   ij     ⁡     (     X   t   ′     )         p   ⁡     (     X   t   ′     )                 {                 U   ⁡     (   t   )       +   1     ≤   i   ≤   U     ,                     V   ⁡     (   t   )       +   1     ≤   j   ≤   V     ,                             (   4   )               
where
 
                          z   t          =         ∑     i   =   1       U   ⁡     (   t   )         ⁢       ∑     j   =   1       V   ⁡     (   t   )         ⁢       z     tij   ,       ⁢           ⁢     P   ⁡     (     X   t     )             =       ∑     i   =   1       U   ⁡     (   t   )         ⁢       ∑     j   =   1       V   ⁡     (   t   )         ⁢       P   ij     ⁡     (     X   t     )               ,           (   5   )               
and X t ′ is the last estimate of X t . Thus, from (4) it may be seen that expected measurements exist for all cells, even those truncated in the observation. After taking the quantization limit,  h   2 →0, the H-PMHT auxiliary functions become
 
                   Q     t   ⁢           ⁢   π       =       ∑     k   =   0     M     ⁢       [       ∑     i   =   1     U     ⁢       ∑     j   =   1     V     ⁢           z   _     tij         P   ij     ⁡     (     X   t   ′     )         ⁢       ∫       C   ij     ⁡     (   t   )         ⁢         G   k     ⁡     (     u   ,     v   ❘     x   tk   ′         )       ⁢           ⁢     ⅆ   u     ⁢     ⅆ   v               ]     ⁢       (     π   ′     )     tk     ⁢   log   ⁢           ⁢     π   tk                 (   6   )               
and
 
                   Q   kx     =         ∑     t   =   1     T     ⁢              z   t            P   ⁡     (     X   t   ′     )         ⁢   log   ⁢           ⁢     P     Ξ   tk           ⁢     ❘     Ξ       t   -   1     ,   k         ⁢       (       x   tk     ❘     x       t   -   1     ,   k         )     +       ∑     t   =   1     T     ⁢       ∑     i   =   1     U     ⁢       ∑     j   =   1     V     ⁢           π   tk   ′     ⁢     z   tij           P   ij     ⁡     (     X   t   ′     )         ⁢       ∫     C   ij     (   t   )         ⁢         G   k     ⁡     (     u   ,     v   ❘     x   tk   ′         )       ⁢   log   ⁢           ⁢       G   k     ⁡     (     u   ,     v   ❘     x   tk         )       ⁢           ⁢     ⅆ   u     ⁢           ⁢       ⅆ   v     .                             (   7   )               
The density P Ξ     tk     |Ξ     t−1,k   (x tk |x t−1,k ) for t=1, 2, . . . T describes the Markov process for the state of target k.
 
   Let the spectral PDF of target k be denoted by S k (v), so that 
                     ∫   ɛ     ⁢         S   k     ⁡     (   v   )       ⁢           ⁢     ⅆ   v         =   1.           (   8   )               
The spectral PDF is equal to the traditional power spectrum normalized so that its integral over ε is one. Because the target spatial and spectral characteristics are independent by assumption, each component G k (u, v|x tk ) of the sample PDF factors:
   G   k ( u,v|x   tk )= g   k ( u|x   tk ) S   k ( v ),  (9) 
where g k (u|x tk ) is the spatial PDF of component k. Independence enables integrals over C ij (t) to be rewritten as products of integrals, so that
 
                     ∫       C   ij     ⁡     (   t   )         ⁢         G   k     ⁡     (     u   ,     v   |     x   tk   ′         )       ⁢           ⁢     ⅆ   u     ⁢           ⁢     ⅆ   v         =       ∫       E   j     ⁡     (   t   )         ⁢         S   k     ⁡     (   v   )       ⁢           ⁢     ⅆ   v     ⁢       ∫       D   i     ⁡     (   t   )         ⁢         g   k     ⁡     (     u   |     x   tk   ′       )       ⁢           ⁢       ⅆ   u     .                     (   10   )               
and, using the mixture (3) and the definition (2),
 
                     P   ij     ⁡     (     X   t   ′     )       =       ∑     k   =   0     M     ⁢           ⁢       π   tk   ′     ⁢       ∫       E   j     ⁡     (   t   )         ⁢         S   k     ⁡     (   v   )       ⁢           ⁢     ⅆ   v     ⁢       ∫       D   i     ⁡     (   t   )         ⁢         g   k     ⁡     (     u   |     x   tk   ′       )       ⁢           ⁢       ⅆ   u     .                         (   11   )               
Substituting (10) into (6) gives
 
                     Q     t   ⁢           ⁢   π       =       ∑     k   =   0     M     ⁢       [       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ⁢       ∫       D   i     ⁡     (   t   )         ⁢         g   k     ⁡     (     u   |     x   tk   ′       )       ⁢           ⁢     ⅆ   u             ]     ⁢     π   tk   ′     ⁢   log   ⁢           ⁢     π   tk           ,           (   12   )               
where
 
                   Ψ   tki     =     (       ∑     j   =   1     V     ⁢           ⁢           z   _     tij     ⁢       ∫       E   j     ⁡     (   t   )         ⁢         S   k     ⁡     (   v   )       ⁢           ⁢     ⅆ   v               P   ij     ⁡     (     X   t   ′     )           )             (   13   )               
is analogous to a normalized matched filter output for target k on spatial cell i at time t, and P ij (X t ′) is given in (11). Similarly, (7) becomes
 
                   Q   kx     =         ∑     t   =   1     T     ⁢           ⁢              Z   t            P   ⁡     (     X   t   ′     )         ⁢   log   ⁢           ⁢       p       Ξ     t   ,   k       |     Ξ       t   -   1     ,   k           ⁡     (       x   tk     |     x       t   -   1     ,   k         )           +       ∑     t   =   1     T     ⁢           ⁢       π   tk   ′     ⁢       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ×       ∫       D   i     ⁡     (   t   )         ⁢         g   k     ⁡     (     u   |     x   tk   ′       )       ⁢           ⁢   log   ⁢           ⁢       g   k     ⁡     (     u   |     x   tk       )       ⁢       ⅆ   u     .                           (   14   )               
There is an additional term in (14), but it is omitted here because it depends on x′ t,k  and not on x t,k , and thus does not influence the M-step of the EM method. It should be noted at this point that it is not necessary to have an analytic expression for S k (v) to utilize (12) and (14). It is sufficient to know the values of the set of integrals
 
             {       ∫       E   j     ⁡     (   t   )         ⁢         S   k     ⁡     (   v   )       ⁢           ⁢     ⅆ   v         }     ,         
j=1, . . . V, for each target k. This vector of spectral cell  38  probabilities is a non-parametric description of the target spectral density sufficient for the problem at hand.
 
   At this stage, specific parametric forms are adopted for the target and measurement processes. For target k, k=1, . . . , M, the process evolution is defined by
 
 p   Ξ     t,k     |Ξ     t−1,k   ( x   tk   |x   t−1,k )= N ( x   tk   ;F   t−1,k   x   t−1,k,   Q   t−1,k )  (15)
 
where N(x; μ, Σ) is the multivariate normal distribution in x with mean μ and covariance Σ. The measurements are characterized by
 
 g   k ( u|x   tk )= N ( u;H   tk   x   tk   ,R   tk ).  (16)
 
The covariance matrix R tk  relates to the spatial extent, or spreading, of the energy about its expected location given by H tk x tk . Estimates of {{circumflex over (π)} tk } are obtained using a Lagrange multiplier technique. The result is
 
                       π   ^     tk     =         π   tk   ′       λ   t       ⁢       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ⁢       ∫       D   i     ⁡     (   t   )         ⁢       N   ⁡     (       u   ;       H   tk     ⁢     x   tk   ′         ,     R   tk       )       ⁢     ⅆ   u                 ,           (   17   )               
where
 
                   λ   t     =         ∑     k   =   0     M     ⁢           ⁢       π   tk   ′     ⁡     [       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ⁢       ∫       D   i     ⁡     (   t   )         ⁢       N   ⁡     (       u   ;       H   tk     ⁢     x   tk   ′         ,     R   tk       )       ⁢           ⁢     ⅆ   u             ]         =       ∑     i   =   1     U     ⁢           ⁢       ∑     j   =   1     V     ⁢           ⁢       z   _     tij                   (   18   )               
The last form follows by taking the sum over k innermost and using Eq. (11).
 
   Estimates for the state variables are obtained by setting the gradient of the auxiliary function Q kX  to zero and solving; however, as in the earlier developments of H-PMHT, an alternative approach is taken because it exploits the Kalman filter as an efficient computational algorithm. The details of the Kalman filter steps are omitted here, however, the synthetic spatial measurements used in the filter for target k now have the form 
                       z   ~     tk     =       1     v   tk       ⁢       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ⁢       ∫       D   i     ⁡     (   t   )         ⁢     u   ⁢           ⁢     N   ⁡     (       u   ;       H   tk     ⁢     x   tk   ′         ,     R   tk       )       ⁢           ⁢     ⅆ   u                 ,           (   19   )               
where
 
                   v   tk     =       ∑     i   =   1     U     ⁢           ⁢       Ψ   tki     ⁢       ∫       D   i     ⁡     (   t   )         ⁢           ⁢       N   ⁡     (       u   ;       H   tk     ⁢     x   tk   ′         ,     R   tk       )       ⁢           ⁢       ⅆ   u     .                     (   20   )               
The synthetic process and measurement noise covariance matrices used in conjunction with this synthetic measurement are respectively given by
 
                       Q   ~     tk     =         P   ⁡     (     X     t   -   1     ′     )              Z     t   +   1              ⁢     Q   tk         ,           ⁢     0   ≤   t   ≤     T   -   1               (   21   )               
and
 
   
     
       
         
           
             
               
                 
                   
                     
                       R 
                       ~ 
                     
                     tk 
                   
                   = 
                   
                     
                       R 
                       tk 
                     
                     
                       
                         π 
                         tk 
                         ′ 
                       
                       ⁢ 
                       
                         v 
                         tk 
                       
                     
                   
                 
                 , 
                 
                     
                 
                 ⁢ 
                 
                   1 
                   ≤ 
                   t 
                   ≤ 
                   T 
                 
               
             
             
               
                 ( 
                 22 
                 ) 
               
             
           
         
       
     
   
   Let {π tk   l } be the set of estimated mixing proportions and {x tk   l } and {R tk   l } define the signal states and width parameters at the l-th EM iteration. For simplicity and robustness, assume that {π tk   l }={π k   l } and {R tk   l }={R k   l } for all t=1, . . . , T in the batch of scans of the sensor space  35 . These restrictions, tantamount to statistical stationarity, are most often reasonable over the data intervals of interest. Further, since the spectral density is never itself required, we will denote the needed integrals by 
   
     
       
         
           
             s 
             kj 
           
           = 
           
             
               ∫ 
               
                 
                   E 
                   j 
                 
                 ⁡ 
                 
                   ( 
                   t 
                   ) 
                 
               
             
             ⁢ 
             
               
                 
                   S 
                   k 
                 
                 ⁡ 
                 
                   ( 
                   v 
                   ) 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   ⅆ 
                   v 
                 
                 . 
               
             
           
         
       
     
   
   The method as described below is illustrated in the flow chart in  FIG. 4 . At the beginning of the method (the 0-th iteration), the mixing proportions {π k   (0) } are initialized so that π k   (0 &gt;0 and π 0   (0) +π 1   (0) + . . . +π M   (0) =1. The signal state sequences x k   (0) ={x 1k   (0) , . . . , x tk   (0) , . . . , x Tk   (0) } are initialized with nominal values for k=1, . . . , M, and the signal widths {R 1   (0) , R 2   (0) , . . . , R M   (0) } are set nominally above the expected signal widths so that the tracks are better able to “see” nearby energy when poorly initialized. The simple case of x k   (0) ={x 0,k   (0) , . . . , x 0,k   (0) , . . . , x 0,k   (0) }, stationary target), has proven an effective starting point in many cases. 
   The process covariance matrices Q t ={Q t,1 , Q t,2 , . . . , Q t,M } are initialized with values tailored to the problem at hand so as to compromise between smooth tracking and the ability to follow through aberrant behavior. Typically it is assumed that the process covariance matrices are constant over time Q t =Q={Q 1 , Q 2 , . . . , Q M }. In order to get the iterative estimator started, initial values are also required for the target state spectral distributions S={S 1 , S 2 , . . . , S M }. The simple case of 
             s   k     =     {       1   v     ,     1   v     ,   ⋯   ⁢           ,     1   v       }           
has proven an effective starting point for estimating the spectra of spatially isolated targets. The above described initialization of target parameters is step  50  in  FIG. 4 .
 
   For iterations l=1, 2, . . . , the following quantities are computed:
         1. Component spatial cell probabilities for t=1, . . . , T, i=1, . . . , U, and k=0, 1, . . . , M:       

   
     
       
         
           
             
               
                 
                   
                     P 
                     ki 
                     
                       ( 
                       l 
                       ) 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       X 
                       t 
                     
                     ) 
                   
                 
                 = 
                 
                   { 
                   
                     
                       
                         
                           
                             1 
                             / 
                             U 
                           
                           , 
                         
                       
                       
                         
                           k 
                           = 
                           0 
                         
                       
                     
                     
                       
                         
                           
                             
                               ∫ 
                               
                                 D 
                                 i 
                               
                             
                             ⁢ 
                             
                               
                                 N 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       τ 
                                       ; 
                                       
                                         
                                           H 
                                           tk 
                                         
                                         ⁢ 
                                         
                                           x 
                                           tk 
                                           
                                             ( 
                                             
                                               l 
                                               - 
                                               1 
                                             
                                             ) 
                                           
                                         
                                       
                                     
                                     , 
                                     
                                       R 
                                       tk 
                                       
                                         l 
                                         - 
                                         1 
                                       
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 ⅆ 
                                 τ 
                               
                             
                           
                           , 
                         
                       
                       
                         
                           k 
                           ≠ 
                           0 
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 23 
                 ) 
               
             
           
         
       
     
       
       
         
           2. Component spatial/spectral cell probabilities for t=1, . . . , T and i=1, . . . , U, j=1, . . . V, and k=0, 1, . . . , M:
 
 P   kij   (l) ( X   t )= P   ki   (l) ( X   t ) S   kj .  (24)
 
           3. Total spatial/spectral cell probabilities for t=1, . . . , T and i=1, . . . , U, j=1, . . . , V: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   
                     P 
                     ij 
                     
                       ( 
                       l 
                       ) 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       X 
                       t 
                     
                     ) 
                   
                 
                 = 
                 
                   
                     ∑ 
                     
                       k 
                       = 
                       0 
                     
                     M 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       π 
                       tk 
                       
                         ( 
                         
                           l 
                           - 
                           1 
                         
                         ) 
                       
                     
                     ⁢ 
                     
                       
                         
                           P 
                           kij 
                           
                             ( 
                             l 
                             ) 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             X 
                             t 
                           
                           ) 
                         
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 25 
                 ) 
               
             
           
         
       
     
       
       
         
           4. Total scan probabilities for t=1, . . . , T: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   P 
                   
                     ( 
                     
                       X 
                       t 
                     
                     ) 
                   
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                     
                     U 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         j 
                         = 
                         1 
                       
                       V 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           P 
                           ij 
                           
                             ( 
                             l 
                             ) 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             X 
                             t 
                           
                           ) 
                         
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 26 
                 ) 
               
             
           
         
       
     
       
       
         
           5. Expected sensor space measurement  z   tij  for t=1, . . . , T i=1, . . . , U, and j=1, . . . , V using equation (4), 
           6. Spatial cell first moments for t=1, . . . , T i=1, . . . , U, and k=1, . . . , M: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   μ 
                   tki 
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     ∫ 
                     
                       D 
                       i 
                     
                   
                   ⁢ 
                   
                     τ 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       N 
                       ( 
                       
                         
                           τ 
                           ; 
                           
                             
                               H 
                               tk 
                             
                             ⁢ 
                             
                               x 
                               tk 
                               
                                 ( 
                                 
                                   l 
                                   - 
                                   1 
                                 
                                 ) 
                               
                             
                           
                         
                         , 
                         
                           R 
                           tk 
                           
                             ( 
                             
                               l 
                               - 
                               1 
                             
                             ) 
                           
                         
                       
                       ⁢ 
                       
                           
                       
                       ) 
                     
                     ⁢ 
                     
                       
                         ⅆ 
                         τ 
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 27 
                 ) 
               
             
           
         
       
     
       
       
         
           7. Relative mode contributions for t=1, . . . , T and k=0, 1, . . . , M: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   v 
                   tk 
                 
                 = 
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                     
                     U 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         j 
                         = 
                         1 
                       
                       V 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           
                             z 
                             tij 
                           
                           ⁢ 
                           
                             
                               P 
                               kij 
                               
                                 ( 
                                 l 
                                 ) 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 X 
                                 t 
                               
                               ) 
                             
                           
                         
                         
                           
                             P 
                             ij 
                             
                               ( 
                               l 
                               ) 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               X 
                               t 
                             
                             ) 
                           
                         
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 28 
                 ) 
               
             
           
         
       
     
       
       
         
           8. Estimated mixing proportions for t=1, . . . , T and k=0, 1, . . . , M: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   π 
                   tk 
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     
                       π 
                       tk 
                       
                         ( 
                         
                           l 
                           - 
                           1 
                         
                         ) 
                       
                     
                     ⁢ 
                     
                       v 
                       tk 
                     
                   
                   
                     
                       ∑ 
                       
                         
                           k 
                           ′ 
                         
                         = 
                         0 
                       
                       M 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         π 
                         
                           tk 
                           ′ 
                         
                         
                           ( 
                           
                             l 
                             - 
                             1 
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         v 
                         tk 
                       
                     
                   
                 
               
             
             
               
                 ( 
                 29 
                 ) 
               
             
           
         
       
     
       
       
         
           9. Synthetic measurements for t=1, . . . , T and k=1, . . . , M: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   
                     z 
                     ~ 
                   
                   tk 
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     1 
                     
                       v 
                       tk 
                     
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         1 
                       
                       U 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         V 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           
                             
                               z 
                               _ 
                             
                             tij 
                           
                           ⁢ 
                           
                             S 
                             kj 
                           
                           ⁢ 
                           
                             μ 
                             tki 
                             
                               ( 
                               l 
                               ) 
                             
                           
                         
                         
                           
                             P 
                             ij 
                             
                               ( 
                               l 
                               ) 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               X 
                               t 
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 30 
                 ) 
               
             
           
         
       
     
       
       
         
           10. Synthetic measurement covariance matrices for t=1, . . . , T and k=1, . . . , M: 
         
       
     
  
   
     
       
         
           
             
               
                 
                   
                     R 
                     ~ 
                   
                   tk 
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     
                       R 
                       tk 
                       
                         ( 
                         
                           l 
                           - 
                           1 
                         
                         ) 
                       
                     
                     
                       
                         π 
                         tk 
                         
                           ( 
                           
                             l 
                             - 
                             1 
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         v 
                         tk 
                       
                     
                   
                   . 
                 
               
             
             
               
                 ( 
                 31 
                 ) 
               
             
           
         
       
     
       
       
         
           11. Synthetic process covariance matrices for t=0, 1, . . . , T−1 and k=1, . . . , M: 
         
       
     
  
                       Q   ~     tk     (   l   )       =           P     (   l   )       ⁡     (     X   t     )              z   t            ⁢     Q   tk         ,           (   32   )               
where Q tk  is treated as a control parameter for the process description, and most commonly Q tk =Q k  for all t=1, . . . , T in the batch.
         12. Estimated spatial states  55  in  FIG. 4  x (l) ={x 01   (l) , . . . , x tk   (l) , . . . , x TM   (l) } for t=0, 1, . . . , T and k=1, . . . , M, using (for computational efficiency) a recursive Kalman smoothing filter, on the synthetic data {tilde over (z)} tk   (l)  with process and measurement matrices corresponding to F tk , {tilde over (Q)} tk   (l) , H tk , {tilde over (R)} tk   (l) .   13. Spatial cell second moments for t=1, . . . , T, i=1, . . . , U. and k=1, . . . , M:       
   
     
       
         
           
             
               
                 
                   σ 
                   tki 
                   
                     ( 
                     l 
                     ) 
                   
                 
                 = 
                 
                   
                     ∫ 
                     
                       D 
                       i 
                     
                   
                   ⁢ 
                   
                     
                       
                         ( 
                         
                           τ 
                           - 
                           
                             
                               H 
                               tk 
                             
                             ⁢ 
                             
                               x 
                               tk 
                               
                                 ( 
                                 
                                   l 
                                   - 
                                   1 
                                 
                                 ) 
                               
                             
                           
                         
                         ) 
                       
                       2 
                     
                     ⁢ 
                     
                       N 
                       ⁡ 
                       
                         ( 
                         
                           
                             τ 
                             ; 
                             
                               
                                 H 
                                 tk 
                               
                               ⁢ 
                               
                                 x 
                                 tk 
                                 
                                   ( 
                                   
                                     l 
                                     - 
                                     1 
                                   
                                   ) 
                                 
                               
                             
                           
                           , 
                           
                             R 
                             tk 
                             
                               ( 
                               
                                 l 
                                 - 
                                 1 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         ⅆ 
                         τ 
                       
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 33 
                 ) 
               
             
           
         
       
     
       
       
         
           14. Average signal width estimates  60  in  FIG. 4  for k=1, . . . , M: 
         
       
     
  
                   R   k     (   l   )       =       (     1       ∑     t   =   1     T     ⁢           ⁢     v   tk         )     ⁢       ∑     t   =   1     T     ⁢           ⁢       ∑     i   =   1     U     ⁢           ⁢       ∑     j   =   1     V     ⁢           ⁢             z   _     tij     ⁢     S   kj     ⁢     σ   tki     (   l   )             P   ij     (   l   )       ⁡     (     X   t     )         .                     (   34   )               
At the completion of iterations l the estimated signal states x (l)  and their width estimates R (l) ={R 1   (l) , R 2   (l) , . . . , R M   (l) } constitute the track estimate output.
         15. Using the track estimate output, compute the average synthetic spectral power  65  in  FIG. 4  for j=1, . . . , V, and k=1, . . . , M:       
   
     
       
         
           
             
               
                 
                   
                     S 
                     ^ 
                   
                   kj 
                 
                 = 
                 
                   
                     ( 
                     
                       1 
                       T 
                     
                     ) 
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         t 
                         = 
                         1 
                       
                       T 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         1 
                         
                           v 
                           tk 
                         
                       
                       ⁢ 
                       
                         
                           ∑ 
                           
                             i 
                             = 
                             1 
                           
                           U 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           
                             
                               
                                 z 
                                 _ 
                               
                               tij 
                             
                             ⁢ 
                             
                               
                                 P 
                                 kij 
                                 
                                   ( 
                                   
                                     l 
                                     + 
                                     1 
                                   
                                   ) 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   X 
                                   t 
                                 
                                 ) 
                               
                             
                           
                           
                             
                               P 
                               ij 
                               
                                 ( 
                                 
                                   l 
                                   + 
                                   1 
                                 
                                 ) 
                               
                             
                             ⁡ 
                             
                               ( 
                               
                                 X 
                                 t 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 35 
                 ) 
               
             
           
         
       
     
       
       
         
           (Note from (35) that 
         
       
     
  
   
     
       
         
           
             
               ∑ 
               
                 j 
                 = 
                 1 
               
               v 
             
             ⁢ 
             
                 
             
             ⁢ 
             
               S 
               kj 
             
           
           ≡ 
           1. 
         
       
     
   
   The resulting combination of processed spatial, signal width and spectral estimates are linked chronologically  70  and displayed as an image on a computer display screen  75  using display methods known in the art. 
   The advantages of the present invention over the prior art are that the resulting method has improved crossing track performance on sources that have some degree of spectral distinction. The present invention also avoids the need for thresholding and peak-picking to produce point measurements. 
   The spectral estimates (35) may be used to initiate this estimator when run on subsequent batches of data. In the preferred embodiment as a new scan is received the oldest scan of the batch is dropped and the estimation method including the steps 1 through 14 (formulae (23) through (35)) as stated above is run in “sliding batch” fashion using the batch length that provides sufficient smoothing without being unnecessarily long. In this case, t in the equations represents the time index within the batch under consideration. The track and spectral estimates from the previous batch are used as initial values to start the iterations as outlined. 
   The target specific spectral estimates (35) constitute outputs unto themselves and can be easily computed for arbitrary track sequences x′ and R′ used in place of x (l)  and R (l) . The resulting spectral estimates have been termed “track conditioned spectral estimates,” and they serve to give a spectral characterization to tracks generated via other means. 
   Obviously many modifications and variations of the present invention may become apparent in light of the above teachings. For example: g k (u/x tk ) may take a parametric form other than the normal density given in (16), g 0 (u) may be other than the uniform density as implied by (23). While it was shown here that the spectrum could be handled in a non-parametric form, the methods are readily extended to treat a parametric spectral description. 
   In light of the above, it is therefore understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described.