Methods for Estimating the Number of Emitters and Their Parameters

1. Field of the Invention

The present invention relates to multiple emitters generally and, more particularly, but not by way of limitation, to novel methods for estimating the number and parameters of emitters.

2. Background Art

The estimation of the number of emitters and their parameters is a topic of considerable interest in the current literature and a number of estimation methods have been developed. See, for example: "Direction of Arrival Estimation Via Extended Phase Interferometry," by Ying-Wah Wu et al., IEEE Transactions on Aerospace and Electronic Systems, Vol. 31, No. 1, January 1995; "Estimation of Amplitude and Phase Parameters of Multicomponent Signals," IEEE Transactions on Signal Processing," by Benjamin Friedlander et al., IEEE Transactions on Signal Processing, Vol. 43, No. 4, April 1995; "Computationally Efficient Angle Estimation for Signals with Known Waveforms," by Jian Li et al., IEEE Transactions on Signal Processing, Vol. 43, No. 9, September 1995; and the references cited therein. These references deal with the processing of a signal impinging on a sensor in order to determine the number and parameters of emitters that are contained in the impinging signal, the result of such processing being the sensor output. The signals themselves may be radio frequency, infrared, optical, acoustic, seismic, or any other type of intermittent or periodic or pulsed emission of energy from one or more emitters to be collected by one or more sensors.

Examples of applications that make use of such sensor outputs typically include the tracking of a multiplicity of emitters by one or more sensors, as in the tracking of a multiplicity of submarines by a field of sonobuoys tuned to detect sonar emissions from the submarines. Algorithms for such tracking are in general use as described in "Multi-Target Tracking with Radar Application," by Samuel Blackman, Artech House, Norwood, Mass., copyright 1986, and several other current texts. Such tracking algorithms have a computational load that is dependent on the number of input vectors presented to the tracking algorithm by the sensors over a unit of time. In dense emitter environments, the computational load required of such tracking algorithms precludes the real time operation of these algorithms.

Accordingly, it is a principal object of the present invention to provide methods for operating on the output of a sensor receiving signals from emitters, which output is collected over several units of time, in order to estimate the number and parameters of emitters in the collected sensor output and reduce the number of input vectors per unit of time that are presented to a tracking algorithm.

It is a further object of the invention to provide such methods that are accurate and efficient.

Other objects of the present invention, as well as particular features, elements, and advantages thereof, will be elucidated in, or be apparent from, the following description and the accompanying drawing figures.

The present invention achieves the above objects, among others, by providing, in a preferred embodiment, a method for estimating the number of emitters of intermittent or periodic pulsed energy waveforms received at a sensor location and the parameters of said intermittent or periodic energy waveforms, comprising: accumulating over a selected period of time output of said sensor, said output consisting of a set of vectors, each of said vectors consisting of at least one element of parameters of direction, frequency, and amplitude, and which element(s) specify(ies) a point in a space defined by said element(s) for each of said waveforms received when at least said parameter of direction is included; dividing said space into a plurality of bins, considering each said vector in turn, and apportioning weight of said each said vector to be accumulated within at least some of said plurality of bins; determining at the end of said selected period of time which of said bins exceed a threshold strength and are, therefore, said emitters; and calculating said at least one element of parameters for each said emitter as appropriate for said vector elements selected for said space.

The file of this patent contains at least one drawing executed in colors.

Understanding of the present invention and the various aspects thereof will be facilitated by reference to the accompanying drawing figures, submitted for purposes of illustration only and not intended to define the scope of the invention, on which:

FIGS. 1-7 illustrate results achieved by the methods of the present invention.

In the problem considered by the method of the present invention, emitters radiate on a single frequency only, the signal shape is known, but the arrival angle, center frequency, amplitude, and time of arrival are all unknown. The signals are of very short duration, and so the simultaneous arrival of signals at the receiver is a rare event. The receiver processes the observed signal assuming that only one emission might be present (using, for example, a matched filer). If an emission is detected, the receiver adds a vector of parameter estimates to a buffer for further processing by the sorting method. At the end of the collection interval, the set of vectors is processed by the sorting method, and the emitter report is passed on for further processing. While in the method of the present invention, the arrival angle is expressed in terms of azimuth, the arrival angle could be expressed in terms of elevation, or the arrival angle could be expressed in terms of both azimuth and elevation.

The development of the method starts with the assumption that the probability density function for the receiver parameter estimates is unimodal and that the estimates are almost independent. Also, it is assumed that the variances of the estimates may be approximated by the variance associated with the minimum detectable signal level. Since the duration of signal collection is short, it is assumed that the amplitude from an emitter will be constant over the collection interval. Under this assumption, amplitude can be used to discriminate between different emitters and different propagation paths from the same emitter. The set of vectors passed from the receiver contains no "noise". There may be vectors in the set that are a result of interference from own ship, multipath from an emitter producing direct path vectors, and "outliers" resulting from interference between near simultaneous arrivals of two emitter signals at the receiver. "Outliers" are rare and tend to drop out after thresholding. any persistent emission, whether from own ship, multipath, or a direct path emitter, will be reported as an emitter for further processing.

Let X.sup.N ={x.sub.1, . . . , x.sub.N } represent the set of output vectors from the receiver. These vectors may be partitioned into disjoint subsets where each subset represents the collection of vectors from an emitter. It is assumed that the probability density function of each subset of vectors is unimodal. The number of emitters, emitter parameters, and partitioning of the set of vectors into emitter subsets is unknown. Hence, the probability density function of any one vector in the set of output vectors is a mixture of unimodal densities where the constants in the mixture and the shapes of each unimodal density are unknown. In the 3D method, a function proportional to the Parzen estimate of the density function for any single output vector is constructed. The local maxima of this function are used to obtain the number of emitters and their associated parameters.

Consider a set of vectors, x.sub.1, . . . x.sub.N where x.sub.n =›.alpha..sub.n .omega..sub.n A.sub.n t.sub.n !.sup.T and where .alpha..sub.n, .omega..sub.n, A.sub.n, and t.sub.n represent azimuth, frequency, amplitude, and time of arrival for the n.sup.th hit. Ignoring arrival time, the remaining entries in each vector specify a point in three dimensions. This three dimensional space is divided up into "bins" whose centers are spaced 1-.sigma. apart in each dimension, where .sigma. is the standard deviation of a parameter represented by a corresponding dimension at a signal-to-noise ratio that is a munimum detectable level (i.e., a signal-to-noise ratio that yields a detection probability of 0.5 at a specified sensor false alarm probability). The closest bin to the n.sup.th hit is found, and all bins within 3-.sigma. of this bin are updated in a fashion similar to computing the Parzen estimate of the multivariate distribution with a Normal kernel. Let {i.sub.n } represent the set of all indices for bins to be updated for the n.sup.th hit. Let i=(i,j,k) where i, j, and k are the indices for azimuth, frequency, and amplitude, respectively. The 3D method computes the function f(.alpha.,.omega.,A) sampled at the bin centers as ##EQU1## for (i,j,k).epsilon..orgate..sup.N {i.sub.n }. f(.alpha.,.omega.,A) is zero for all other (i,j, k). The Normal kernel in (1) is not normalized by the variance, as as in the Parzen estimate of the multivariate probability density function. Note that a vector that is centered in a bin contributes "one" to the result for the function and a relatively small amount to the surrounding bins. After the set of vectors is processed, the sampled values of f can be interpreted as a weighted average of the number of hits in the vicinity of each bin. Those samples that lie above a threshold are processed further. Since it is assumed that the probability density functions for the measurements from each emitter are all unimodal, the method tests those bins above threshold to determine if they are local maxima. The estimates of the emitter parameters are obtained by a parabolic interpolation about the local maxima. Let (i.sub.l,j.sub.l,k.sub.l),l=1, . . . , L represent the indices of the bins that are greater than the threshold, .eta.. Let {m.sub.l } represent the set of indices within .kappa..sigma. cube centered at the l.sup.th bin, and excluding the l.sup.th bin. Then the set of bins that contains emitters is given by ##EQU2## where f.sub.l represents f evaluated at the l.sup.th bin. Once the emitter bins have been found, the parameters are obtained using a parabolic interpolation. For example, to estimate azimuth, let f.sub.p-1 =f((i.sub.p -1).sigma..sub..alpha.,j.sub.p .sigma..sub..omega.,k.sub.p .sigma..sub.A), f.sub.p =f(i.sub.p .sigma..sub..alpha.,j.sub.p .sigma..sub..omega.,k.sub.p .sigma..sub.A), and f.sub.p+1 =f((i.sub.p +1).sigma..sub..alpha.,j.sub.p .sigma..sub..omega.,k.sub.p .sigma..sub.A). Then the estimate of the p.sup.th emitter azimuth is given by ##EQU3##

Similar results apply for the estimates of frequency and amplitude. Since the function f is proportional to the probability density function for the p.sup.th emitter in the neighborhood of the p.sup.th bin (ignoring the distortion from the other pdf's), the estimates can be interpreted as maximum a posteriori estimates of the parameters for the p.sup.th emitter.

The method occasionally produces two identical emitters, due to the use of equality in equation (2). Without equality in equation (2), the method would fail to report some emitters. The constants used in the description of the method (e.g., .sigma..sub..alpha., .sigma..sub..omega., .sigma..sub.A, .eta., and .kappa.) are parameters that can be tuned to optimize performance. Also, the exponential in (1) can be replaced with a suitable positive function to speed up the computation or optimize performance. The algorithm presented here will be referred to as the 3D method. In the 2D method described below, a slight modification to this algorithm is presented to speed up the processing required for very large data sets. This modification results in a slight degradation in performance, but a significant increase in speed.

In this method, one dimension contains a very small number of bins and each bin in this dimension is treated like a plane in testing for the local maximum. There are three bins (planes) in amplitude for this method, the closest bin for the n.sup.th hit is found in each amplitude plane, and all bins within 3-.sigma. in azimuth and frequency are updated. Let {i.sub.n,k } represent the set of all indices for bins to be updated for the n.sup.th hit in the k.sup.th plane. The 2D method computes the function ##EQU4## for k=1, 2, 3. Each plane is processed separately, and the emitters for each plane are merged to form the final set of reported emitters. The function g.sub.k is given by ##EQU5## Where A.sub.max is the maximum of the amplitude, .left brkt-top.x.right brkt-top. is the smallest integer greater than x, and .left brkt-bot.x.right brkt-bot. is the greatest integer less than x. Let (i.sub.l,j.sub.l, k), l=1, . . . L.sub.k represent the indices of the L.sub.k bins that are greater than the threshold, .eta., in the k.sup.th amplitude plane. Let {m.sub.l,k } represent the set of indices within a .kappa..sigma. square centered at the l.sup.th bin in the k.sup.th amplitude plane, excluding the l.sup.th bin. Then the set of bins in the k.sup.th plane that contain emitters is ##EQU6## where f.sub.l,k represents f evaluated at the l.sup.th bin in the k.sup.th amplitude plane. Once the emitter bins have been found, azimuth and frequency are estimated using parabolic interpolation, as in equation (3). To obtain the final set of emitters, redundant emitters in each plane are removed first. Next, emitters in adjacent planes are compared using the Mahalinobis distance ##EQU7## where k.sub.1 and k.sub.2 are adjacent amplitude planes and every estimate in the k.sub.1 plane is compared with every estimate in the k.sub.2 plane. When estimates in the adjacent planes come from the same emitter, d.sup.2.sub.k.sbsb.1.sub.k.sbsb.2 is central chi-squared with two degrees of freedom (assuming Gaussian errors on the estimates). If d.sup.2.sub.k.sbsb.1.sub.k.sbsb.2 <9.21 (a miss probability of 0.01), then the emitter with the smaller f is dropped. Note that emitters in planes 1 and 3 that satisfy this condition are unaffected. This permits the detection of a distant emitter at nearly the same azimuth and frequency as a nearby emitter.

Both methods were evaluated using over 100 data sets representing typical target densities. The data varies from sparse, containing about 50 hits, to dense, containing over 1000 hits. Azimuth is an integer in the range -128 to +128, frequency is an integer in the range -30 to +30, and amplitude is an integer in the range 0 to A.sub.max =256. The parameters used for the 2D method were .sigma..sub..alpha. =2, .sigma..sub..omega. =3, .kappa.=2.0, and .eta.=2.0. An azimuth of zero corresponds to the boresight of the receive antenna and a frequency of zero corresponds to the center of the band to which the receiver is tuned.

There were only minor variations between the emitter reports for the 2D and 3D methods in sparse environments. A typical example of the 2D method in a sparse environment is presented on FIGS. 1-3, which shows the amplitude planes, those hits that have amplitudes that lie in the plane, and emitters reported in each plane.

The emitter report for this example is present in Table I.

An example of the 2D method in a dense environment is shown on FIGS. 4-6. The threshold setting for this example was .eta.=4.0. The emitter report for this example is presented in Table II.

The dense data set was also evaluated using the 3D method. The results are presented on FIG. 7 and in Table III. To obtain FIG. 7, all of the azimuth and frequency bins were added in the plane of the peak amplitude response for each emitter.

Comparison of Tables II and III reveals that the emitter at an azimuth and frequency of -33.7491 and -17.8938 in the 2D method has been resolved into two emitters, one at (-33.8889, -18.3197) with an amplitude of 56.2207 and another at (-33.4136, -17.4107) with an amplitude of 93.5981. The remaining emitters are close to the 2D emitter report.

The data show that the 2D method emitter reports are very close to the reports of the 3D method for typical target densities. There is, in fact, only one data set of the 100 in the test where the 2D method appears to miss an emitter found by the 3D method. It is expected that the 2D method can process the receiver buffer in real time, even in dense emitter environments.

It will thus be seen that the objects set forth above, among those elucidated in, or made apparent from, the preceding description, are efficiently attained and, since certain changes may be made in the above methods without departing from the scope of the invention, it is intended that all matter contained in the above description or shown on the accompanying drawing figures shall be interpreted as illustrative only and not in a limiting sense.

It is also to be understood that the following claims are intended to cover all of the generic and specific features of the invention herein described and all statements of the scope of the invention which, as a matter of language, might be said to fall therebetween.