Method and system for tracking multiple regional objects by multi-dimensional relaxation

[0002] This application is a continuation of U.S. patent application Ser. No. 08/682,904 filed Jul. 16, 1996, now U.S. Pat. No. 5,959,574 which is a continuation-in-part of U.S. patent application Ser. No. 08/404,024, filed Mar. 14, 1995 and issued Jul. 16, 1996 as U.S. Pat. No. 5,537,119, which is a continuation-in-part of U.S. patent application Ser. No. 08/171,327 filed Dec. 21, 1993, now U.S. Pat. No. 5,406,289.

[0003] a. Field of the Invention

[0004] The invention relates generally to computerized techniques for processing data obtained from radar to track multiple discrete objects.

[0005] b. Description of the Background

[0006] There are many situations where the courses of multiple objects in a region must be tracked. Typically, radar is used to scan the region and generate discrete images or “snapshots” based on sets of returns or observations. In some types of tracking systems, all the returns from any one object are represented in an image as a single point unrelated to the shape or size of the objects. “Tracking” is the process of identifying a sequence of points from a respective sequence of the images that represents the motion of an object. The tracking problem is difficult when there are multiple closely spaced objects because the objects can change speed and direction rapidly and move into and out of the line of sight for other objects. The problem is exacerbated because each set of returns may result from noise as well as echoes from the actual objects. The returns resulting from the noise are also called false positives. Likewise, the radar will not detect all echoes from the actual objects and this phenomena is called a false negative or “missed detect” error. For tracking airborne objects, a large distance between the radar and the objects diminishes the signal to noise ratio so the number of false positives and false negatives can be high. For robotic applications, the power of the radar is low and as a result, the signal to noise ratio can also be low and the number of false positives and false negatives high.

[0007] In view of the proximity of the objects to one another, varied motion of the objects and false positives and false negatives, multiple sequential images should be analyzed collectively to obtain enough information to properly assign the points to the proper tracks. Naturally, the larger the number of images that are analyzed, the greater the amount of information that must be processed.

[0008] While identifying the track of an object, a kinematic model may be generated describing the location, velocity and acceleration of the object. Such a model provides the means by which the future motion of the object can be predicted. Based upon such a prediction, appropriate action may be initiated. For example, in a military application there is a need to track multiple enemy aircraft or missiles in a region to predict their objective, plan responses and intercept them. Alternatively, in a commercial air traffic control application there is a need to track multiple commercial aircraft around an airport to predict their future courses and avoid collision. Further, these and other applications, such as robotic applications, may use radar, sonar, infrared or other object detecting radiation bandwidths for tracking objects. In particular, in robotic applications reflected radiation can be used to track a single object which moves relative to the robot (or vice versa) so the robot can work on the object.

[0009] Consider the very simple example of two objects being tracked and no false positives or false negatives. The radar, after scanning at time t₁, reports objects at two locations in a first observation set. That is, the radar returns a set of two observations {o₁₁, o₁₂}. At time t₂the radar returns a similar set of two observations {o₂₁, o₂₂) from a second observation set. Suppose from prior processing that track data for two tracks T₁and T₂includes the locations at t₀of two objects. Track T₁may be extended through the points in the two sets of observations in any of four ways, as may track T₂. The possible extensions of T₁can be described as: {T₁, o₁₁, o₂₁}, {T₁, o₁₁, o₂₂}, {T₁, o₁₂, o₂₁} and {T₁, o₁₂, o₂₂}. Tracks can likewise be extended from T₂in four possible ways, including {T₂, o₁₂, o₂₁}. FIG. 1 illustrates these five (out of eight) possible tracks (to simplify the problem for purposes of explanation). The five track extensions are labeled h₁₁, h₁₂, h₁₃, h₁₄, and h₂₁wherein h₁₁is derived from {T₁, o₁₁, o₂₁}, h₁₂is derived from {T₁, o₁₁, o₂₂}, h₁₃is derived from {T₁, o₁₂, o₂₁}, h₁₄is derived from {T₁, o₁₂, o₂₂}, and h₂₁is derived from {T₂, o₁₁, o₂₁}. The problem of determining which such track extensions are the most likely or optimal is hereinafter known as the assignment problem.

[0010] It is known from prior art to determine a figure of merit or cost for assigning each of the points in the images to a track. The figure of merit or cost is based on the likelihood that the point is actually part of the track. For example, the figure of merit or cost may be based on the distance from the point to an extrapolation of the track. FIG. 1 illustrates costs δ₂₁and δ₂₂for hypothetical extension h₂₁and modeled target characteristics. The function to calculate the cost will normally incorporate detailed characteristics of the sensor, such as probability of measurement error, and track characteristics, such as likelihood of track maneuver.

[0011] FIG. 2 illustrates a two by two by two matrix, c, that contains the costs for each point in relation to each possible track. The cost matrix is indexed along one axis by the track number, along another axis by the image number and along the third axis by a point number. Thus, each position in the cost matrix lists the cost for a unique combination of points and a track, one point from each image. FIG. 2 also illustrates a {0, 1} assignment matrix, z, which is defined with the same dimensions as the cost matrix. Setting a position in the assignment matrix to “one” means that the equivalent position in the cost matrix is selected into the solution. The illustrated solution matrix selects the {h₁₄, h₂₁} solution previously described. Note that for the above example of two tracks and two snapshots, the resulting cost and assignment matrices are three-dimensional. As used in this patent application, the term “dimension” means the number of axes in the cost or assignment matrix while size refers to the number of elements along a typical axis. The costs and assignments have been grouped in matrices to facilitate computation.

[0012] A solution to the assignment problem satisfies two constraints—first, the sum of the associated costs for assigning points to a track extension is minimized and, second, if no false positives or false negatives exist, then each point is assigned to one and only one track.

[0013] When false positives exist, however, additional hypothetical track extensions incorporating the false positives will be generated. Further note that the random locations of false positives will, in general, not fit well with true data and such additional hypothetical track extensions will result in higher costs. Also note that when false negative errors exist, then the size of the cost matrix must grow to include hypothetical track extensions formulated with “gaps” (i.e., data omissions where there should be legitimate observation data) for the false negatives. Thus, the second criteria must be weakened to reflect false positives not being as signed and also to permit the gap filler to be multiply assigned. With hypothetical cost calculated in this manner then the foregoing criteria for minimization will to materialize the false negatives and avoid the false positives.

[0014] For a three-dimensional problem, as is illustrated in FIG. 1, but with N₁(initial) tracks, N₂observations in scan 1, N₃observations in scan 2, false positives and negatives assumed, the assignment problem can be formulated as: $\begin{array}{cc}\begin{array}{cccc}\left(a\right)& \text{Minimize:}& \sum _{i_1=0}^{N_1}\sum _{i_2=0}^{N_2}\sum _{i_3=0}^{N_3}{C}_{i_1{i}_2{i}_3}{z}_{i_1{i}_2{i}_3}& \\ \left(b\right)& \text{Subject to:}& \sum _{i_2=1}^{N_2}\sum _{i_3=1}^{N_3}{z}_{i_1{i}_2{i}_3}=1,& i_1=1,\dots N_1,\\ \left(c\right)& & \sum _{i_1=1}^{N_1}\sum _{i_3=1}^{N_3}{z}_{i_1{i}_2{i}_3}\le 1,& i_2=1,\dots N_2,\\ \left(d\right)& & \sum _{i_1=1}^{N_1}\sum _{i_2=1}^{N_2}{z}_{i_1{i}_2{i}_3}\le 1,& i_3=1,\dots N_3,\\ \left(e\right)& & z_{i_1{i}_2{i}_3}\in \left\{0,1\right\}& \forall z_{i_1{i}_2{i}_3},\end{array}& \left(0.1\right)\end{array}$

[0015] where “c” is the cost and “z” is a point or observation assignment, as in FIG. 2.

[0016] The minimization equation or equivalently objective function (0.1)(a) specifies the sum of the element by element product of the c and z matrices. The summation includes hypothesis representations z_iiiwith observation number zero being the gap filler observation. Equation (0.1)(b) requires that each of the tracks T₁, . . . , T_Nbe extended by one and only one hypothesis. Equation (0.1)(c) relates to each point or observation in the first observation set and requires that each such observation, except the gap filler, can only associate with one track but, because of the “less than” condition, it might not associate with any track. Equation (0.1)(d) is like (0.1)(c) except that it is applicable to the second observation set. Equation (0.1)(e) requires that elements of the solution matrix z be limited to the zero and one values.

[0017] The only known method to solve Problem Formulation (0.1) exactly is a method called “Branch and Bound.” This method provides a systematic ordering of the potential solutions so that solutions with a same partial solution are accessible via a branch of a tree describing all possible solutions whereby the cost of unexamined solutions on a branch are incrementally developed as the cost for other solutions on the branch are determined. When the developing cost grows to exceed the previously known minimal cost (i.e., the bound) then enumeration of the tree branch terminates. Evaluation continues with a new branch. If evaluation of the cost of a particular branch completes, then that branch has lower cost than the previous bound so the new cost replaces the old bound. When all possible branches are evaluated or eliminated then the branch that had resulted in the last used bound is the solution. If we assume that N₁=N₂=N₃=n and that branches typically evaluate to half there full length, then workload associated with “branch and bound” is proportional to (n!|n/2!)². This workload is unsuited to real time evaluation.

[0018] The Branch and Bound algorithm has been used in past research on the Traveling Salesman Problem. Messrs. Held and Karp showed that if the set of constraints was relaxed by a method of Lagrangian Multipliers (described in more detail below) then tight lower bounds could be developed in advance of enumerating any branch of the potential solution. By initiating the branch and bound algorithm with such a tight lower bound, significant performance improvements result in that branches will typically evaluate to less than half their full length.

[0019] Messrs. Frieze and Yadagar in dealing with a problem related to scheduling combinations of resources, as in job, worker and work site, showed that Problem Formulation (0.1) applied. They further described a solution method based upon an extension of the Lagrangian Relaxation previously mentioned. The two critical extensions provided by Messrs. Frieze and Yadagar were: (1) an iterative procedure that permitted the lower bound on the solution to be improved (by “hill climbing” described below) and (2) the recognition that when the lower bound of the relaxed problem was maximized, then there existed a method to recover the solution of the non-relaxed problem in most cases using parameters determined at the maximum. The procedures attributed to Messrs. Frieze and Yadagar are only applicable to the three-dimensional problem posed by Problem Formulation (0.1) and where the cost matrix is fully populated. However, tracking multiple airborne objects usually requires solution of a much higher dimensional problem.

[0020] FIGS. 1 and 2 illustrate an example where “look ahead” data from the second image improved the assignment accuracy for the first image. Without the look ahead, and based only upon a simple nearest neighbor approach, the assignments in the first set would have been reversed. Problem Formulation (0.1) and the prior art only permit looking ahead one image. In the prior art it was known that the accuracy of assignments will improve if the process looks further ahead, however no practical method to optimally incorporate look ahead data existed. Many real radar tracking problems involve hundreds of tracks, thousands of observations per observation set and matrices with dimensions in the range of 3 to 25 including many images of look ahead.

[0021] It was also known that the data assignment problem may be simplified (without reducing the dimension of the assignment problem) by eliminating from consideration for each track those points which, after considering estimated limits of speed and turning ability of the objects, could not physically be part of the track. One such technique, denoted hereinafter the “cone method,” defines a cone as a continuation of each previously determined track with the apex of the cone at the end of the previously defined track. The length of the cone is based on the estimated maximum speed of the object and the size of the arc of the cone is based on the estimated maximum turning ability of the object. Thus, the cone defines a region outside of which no point could physically be part of the respective track. For any such points outside of the cones, an infinite number could be put in the cost matrix and a zero could be preassigned in the assignment matrix. It was known for the tracking problem that these elements will be very common in the cost and selection matrices (so these matrices are “sparse”).

[0022] It was also known in the prior art that one or more tracks which are substantially separated geographically from the other tracks can be separated also in the assignment problem. This is done by examining the distances from each point to the various possible tracks. If the distances from one set of points are reasonably short only in relation to one track, then they are assigned to that track and not further considered with the remainder of the points. Similarly, if a larger group of points can only be assigned to a few tracks, then the group is considered a different assignment problem. Because the complexity of assignment problems increases dramatically with the number of possible tracks and the total number of points in each matrix, this partitioning of the group of points into a separate assignment problem and removal of these points from the matrices for the remaining points, substantially reduces the complexity of the overall assignment problem.

[0023] A previously known Multiple Hypothesis Testing (MHT) algorithm (see Chapter 10 of S. S. Blackman. Multiple-Target Tracking with Radar Applications. Artech House, Norwood, Mass., 1986.) related to formulation and scoring. The MHT procedure describes how to formulate the sparse set of all reasonable extension hypothesis (for FIG. 1 the set {h₁₁. . . h₂₄}) and how to calculate a cost of the hypothesis {T_i, o_1j, o^2k} based upon the previously calculated cost for hypothesis (T_i, o_1j}. The experience with the MHT algorithm, known in the prior art, is the basis for the assertion that look ahead through k sets of observations results in improved assignment of observations from the first set to the track.

[0024] In theory, the MHT procedure uses the extendable costing procedure to defer assignment decision until the accumulated evidence supporting the assignment becomes overwhelming. When it makes the assignment decision, it then eliminates all potential assignments invalidated by the decision in a process called “pruning the tree.” In practice, the MHT hypothesis tree is limited to a fixed number of generations and the overwhelming evidence rule is replaced by a most likely and feasible rule. This rule considers each track independently of others and is therefore a local decision rule.

[0025] A general object of the present invention is to provide an efficient and accurate process for assigning each point object in a region from multiple images to a proper track and then taking an action based upon the assignments.

[0026] A more specific object of the present invention is to provide a technique of the foregoing type which determines the solution of a k-dimensional assignment problem where “k” is greater than or equal to three.

[0027] The present invention relates to a method and apparatus for tracking objects. In particular, the present invention tracks movement or trajectories of objects by analyzing radiation reflected from the objects, the invention being especially useful for real-time tracking in noisy environments.

[0028] In providing such a tracking capability, a region containing the objects is repeatedly scanned to generate a multiplicity of sequential images or data observation sets of the region. One or more points (or equivalently observations), in each of the images or observation sets are detected wherein each such observation either corresponds to an actual location of an object or is an erroneous data point due to noise. Subsequently, for each observation detected, figures of merit or costs are determined for assigning the observation to each of a plurality of previously determined tracks. *Afterwards, a first optimization problem is specified which includes:

[0029] (a) a first objective function for relating the above mentioned costs to potential track extensions through the detected observations (or simply observations); and

[0030] (b) a first collection of constraint sets wherein each constraint set includes constraints to be satisfied by the observations in a particular scan to which the constraint set is related. In general, there is a constraint for each observation of the scan, wherein the constraint indicates the number of track extensions to which the observation may belong.

[0031] In particular, the first optimization problem is formulated, generated or defined as an M-dimensional assignment problem wherein there are M constraint sets in the first collection of constraint sets (i.e., there are M scans being examined) and the first objective function minimizes a total cost for assigning observations to various track extensions wherein terms are included in the cost, such that the terms have the figures of merit or costs for hypothesized combinations of assignments of the observations to the tracks. Subsequently, the formulated M-dimensional assignment problem is solved by reducing the complexity of the problem by generating one or more optimization problems each having a lower dimension and then solving each lower dimension optimization problem. That is, the M-dimensional assignment problem is solved by solving a plurality of optimization problems each having a lower number of constraint sets.

[0032] The reduction of the M-dimensional assignment problem to a lower dimensioned problem is accomplished by relaxing the constraints on the points of one or more scans thereby permitting these points to be assigned to more than one track extension. In relaxing the constraints, terms having penalty factors are added into the objective function thereby increasing the total cost of an assignment when one or more points are assigned to more than one track. Thus, the reduction in complexity by this relaxation process is iteratively repeated until a sufficiently low dimension is attained such that the lower dimensional problem may be solved directly by known techniques.

[0033] In one embodiment of the invention, each k-dimensional assignment problem (2<k<M) is iteratively reduced to a (k−1)-dimensional problem until a two-dimensional problem is specified or formulated. Subsequently, the two-dimensional problem formulated is solved directly and a “recovery” technique is used to iteratively recover an optimal or near-optimal solution to each k-dimensional problem from a derived (k−1)-dimensional problem k=2,3,4, . . . M.

[0034] In performing each recovery step (of obtaining a solution to a k-dimensional problem using a solution to a (k−1)-dimensional problem), an auxiliary function is utilized. In particular, to recover an optimal or near-optimal solution to a k-dimensional problem, an auxiliary function, Ψ_k−1(k=4,5, . . . ,M), is specified and a region or domain is determined wherein this function is maximized, whereby values of the region determine the penalty factors of the (k−1)-dimensional problem such that another two-dimensional problem can be formulated which determines a solution to the k-dimensional problem using the penalty factors of the (k−1)-dimensional problem.

[0035] Each Auxiliary function Ψ_k, k=3, . . . ,M−1, is a function of both lower dimensional problems, penalty factors, and a solution at the dimension k at which the penalized cost function is solved directly (typically a two-dimensional problem). Further, in determining, for auxiliary function Ψ_k, a peak region, a gradient of the auxiliary function is determined, and utilized to identify the peak region. Thus, gradients are used for each of the approximation functions Ψ₃,Ψ₄, . . . ,Ψ_M−1which are sequentially formulated and used in determining penalty factors until M−1 is used in determining the penalty factors for the (M−1)-dimensional problem. Subsequently, once the M-dimensional problem is solved (using a two-dimensional problem to go from an (M−1)-dimensional solution to an M-dimensional solution), one or more of the following actions are taken based on the track assignments: sending a warning to aircraft or a ground or sea facility, controlling air traffic, controlling anti-aircraft or anti-missile equipment, taking evasive action, working on one of the objects.

[0036] According to one feature of this first embodiment of the present invention, the following steps are also performed before the step of defining the auxiliary function. A preliminary auxiliary function is defined for each of the lower dimensional problems having a dimension equal or one greater than the dimension at which the penalized cost function is solved directly. The preliminary auxiliary function is a function of lower order penalty values and a solution at the dimension at which the penalized cost function was solved directly. In determining a gradient of the preliminary auxiliary function, step in the direction of the gradient to identify a peak region of the preliminary auxiliary function and determine penalty factors at the peak region. Iteratively repeat the defining, gradient determining, stepping and peak determining steps to define auxiliary functions at successively higher dimensions until the auxiliary function dimension (k−1) is determined. In an alternative second embodiment of the present invention, instead of reducing the dimentionality of the M-dimensional assignment problem by a single dimension at a time, a plurality of dimensions are relaxed simultaneously. This new strategy has the advantage that when the M-dimensional problem is relaxed directly to a two-dimensional assignment problem, then all computations may be performed precisely without utilizing an auxiliary function such as Ψ_kas in the first embodiment. More particularly, the second embodiment solves the first optimization problem (i.e., the M-dimensional assignment problem) by specifying (i.e., creating, generating, formulating and/or defining) a second optimization problem. The second optimization problem includes a second objective function and a second collection of constraint sets wherein:

[0037] a) the second objective function is a combination of the first objective function and penalty factors or terms determined for incorporating (M−m)-constraint sets of the first optimization problem into the second objective function;

[0038] b) the constraint sets of the second collection include only m of the constraint sets of the first collection of constraints,

[0039] wherein 2<m<M−1. Note that, once the second optimization problem has been specified or formulated, an optimal or near-optimal solution is determined and that solution is used in specifying (i.e, creating, generating, formulating and/or defining) a third optimization problem of (M−m) dimensions (or equivalently constraint sets). The third optimization problem is subsequently solved by decomposing it using the same procedure of this second embodiment as was used to decompose the first optimization problem above. Thus, a plurality of instantiations of the third optimization problem are specified, each successive instantiation having a lower number of dimensions, until an instance of the third optimization problem is a two-dimensional assignment problem which can be solved directly. Subsequently, whenever an instance of the third optimization problem is solved, the solution is used to recover a solution to the instance of the first optimization problem from which this instance of the third optimization was derived. Thus, an optimal or near-optimal solution to the original first optimization problem may be recovered through iteration of the above steps.

[0040] As mentioned above, the second embodiment of the present invention is especially advantageous when m=2, since in this case all computations may be performed precisely and without utilizing auxiliary functions.

[0041] Other features and benefits of the present invention will become apparent from the detailed description with the accompanying drawings contained hereinafter.

[0042] FIG. 1 is a graph of images or data sets generated by a scan of a region and possible tracks within the images or data sets according to the prior art.

[0043] FIG. 2 illustrates cost and assignment matrices for the data sets of FIG. 1 according to the prior art.

[0044] FIG. 3 is a block diagram of the present invention.

[0045] FIG. 4 is a flow chart of a process according to the prior art for solving a three-dimensional assignment problem.

[0046] FIG. 5 is a flow chart of a process according to the present invention for solving a k-dimensional assignment problem where “k” is greater than or equal to 3.

[0047] FIG. 6 is a graph of various functions used to explain the present invention.

[0048] FIG. 7 is another block diagram of the present invention for solving the k-dimensional assignment problem where “k” is greater than or equal to three (3).

[0049] FIG. 8 is a flowchart describing the procedure for solving an n-dimensional assignment problem according to the second embodiment of the invention.

[0050] Referring now to the other figures in detail wherein like reference numerals indicate like elements throughout the several views, FIG. 3 illustrates a system generally designated 100 for implementing the present invention. System 100 comprises, for example, a radar station 102 (note sonar, microwave, infrared and other radiation bandwidths are also contemplated) for scanning a region which may be, for example, an aerial region (in aerial surveillance applications) or a work region (in robotic applications) and generating signals indicating locations of objects within the region. The signals are input to a converter 104 which converts the signals to data points or observations in which each object (or false positive) is represented by a single point. The output of the converter is input to and readable by a computer 106. As described in more detail below, the computer 106 assigns the points to respective tracks, and then displays the tracks and extrapolations of the tracks on a monitor 110. Also, the computer 106 determines an appropriate action to take based on the tracks and track extensions. For example, in a commercial application at an airport, the computer can determine if two aircraft being tracked are on a collision course and if so, signal a transmitter 112 to warn each aircraft, or if a scheduled take-off will pose the risk of collision, delay the take-off. For a military application on a ship or base, the computer can determine subsequent coordinates of enemy aircraft and send the coordinates to an antiaircraft gun or missile 120 via a communication channel 122. In a robotic application, the computer controls the robot to work on the proper object or portion of the object.

[0051] The invention generates k-dimensional matrices where k is the number of images or sets of observation data in the look ahead window plus one. Then, the invention formulates a k-dimensional assignment problem as in (0.1) above.

[0052] The k-dimensional assignment problem is subsequently relaxed to a (k−1)-dimensional problem by incorporating one set of constraints into the objective function using a Lagrangian relaxation of this set. Given a solution of the (k−1)-dimensional problem, a feasible solution of the k-dimensional problem is then reconstructed. The (k−1)-dimensional problem is solved in a similar manner, and the process is repeated until it reaches the two-dimensional problem that can be solved exactly. The ideas behind the Lagrangian relaxation scheme are outlined next.

[0053] Consider the integer programming problem $\begin{array}{cc}\begin{array}{cc}\mathrm{Minimize}& v\left(z\right)=c^{T}z\\ \text{Subject to:}& \mathrm{Az}\le b,\\ & \mathrm{Bz}\le d,\\ & z_i\mathrm{is}\mathrm{an}\mathrm{integer}\mathrm{for}i\in I,\end{array}& \left(0.2\right)\end{array}$

[0054] where the partitioning of the constraints is natural in some sense. Given a multiplier vector u≧0, the Lagrangian relaxation of (0.2) relative to the constraints Bz≦d is defined to be $\begin{array}{cc}\begin{array}{cc}\Phi \left(u\right)=\mathrm{Minimize}& \left\{c^{T}z+u^T\left(\mathrm{Bz}-d\right)\right\}\\ \text{Subject to:}& \mathrm{Az}\le b,\\ & z_i\mathrm{is}\mathrm{an}\mathrm{integer}\mathrm{for}i\in I.\end{array}& \left(0.3\right)\end{array}$

[0055] If the constraint set Bz≦d is replaced by Bz=d, the non-negativity constraint on u is removed. =c^Tz+u^T(Bz−d) is the Lagrangian relative to the constraints Bz≦d, and hence the name Lagrangian relaxation. The following fact gives the relationship between the objective functions of the original and relaxed problems.

[0056] FACT A.1. If {overscore (z)} is an optimal solution to (0.2), then Φ(u)≦v({overscore (z)}) for all u≧0. If an optimal solution {circumflex over (z)} of (0.3) is feasible for (0.3), then {circumflex over (z)} an optimal solution for (0.2) and Φ(u)=v({overscore (z)}).

[0057] This leads to the following algorithm:

[0058] Algorithm. Construct a sequence of multipliers {u_k}_k=0^∞ converging to the solution {overscore (u)} of Maximize {Φ(u):u≧0} and a corresponding sequence of feasible solutions {z_k}_k=0^∞ of (0.2) as follows:

[0059] 1. Generate an initial approximation u₀.

[0060] 2. Given u_k, choose a search direction s_kand a search distance α_kso that Φ(u_k+α_ks_k)>Φ(u_k). Update the multiplier estimate u_kby u_k+1=u_k+α_ks_k.

[0061] 3. Given u_k+1and a feasible solution {circumflex over (z)}_k+1(u_k+1) of (3.0), recover a feasible solution z_k+1(u_k+1) of the integer programming problem (0.2).

[0062] 4. Check the termination criteria. If the algorithm is not finished, set k=k+1 and return to Step 2. Otherwise, terminate the algorithm.

[0063] If {overscore (z)} is an optimal solution of (0.2), then

[0064] Φ(u_k)≦Φ({overscore (u)})≦v({overscore (z)})≦v(z_k).

[0065] Since the optimal solution {circumflex over (z)}={circumflex over (z)}(u) of (0.3) is usually not a feasible solution of (0.2) for any choice of the multipliers, Φ({overscore (u)}) is usually strictly less than v({overscore (z)}). The difference v({overscore (z)})−Φ({overscore (u)}) is called the “duality gap,” which can be estimated by comparing the best relaxed and recovered solutions computed in the course of maximizing Φ(u).

[0066] Thus, to eliminate the “<” in the constraint formulations (e.g., (0.1)(a)-(0.1)(d)) that resulted from false positives, the constraint sets have been modified, as indicated above, to incorporate a gap filler in each observation set to account for false positives. Thus, a zero for an index i_p(1≦p≦k) in z_{i. . . i}^kin problem (0.4) below indicates that the track extension represented by z_{i. . . i}^kincludes a false positive in the p^thobservation set. Note that this implies that a hypothesis be formed incorporating an observation with k−1 gap fillers, e.g., Z_{0 . . . 0i0 . . . 0},i_p≠0. Thus, the resulting generalization of Problem Formulation (0.1) without the “less than” complication within the constraints is the following k-Dimensional Assignment Problems in which k≧3: $\begin{array}{cc}\begin{array}{cc}\text{Minimize}& \sum _{i_1=0}^{N_1}\cdots \sum _{i_k=0}^{N_k}{c}_{i_{1\cdots }{i}_k}{z}_{i_{1\cdots }{i}_k}\\ \text{Subject to:}& \sum _{i_2=0}^{N_2}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}=1,i_1=1,\dots N_1,\\ & \sum _{i_1=0}^{N_1}\cdots \sum _{i_{j-1}=0}^{N_{j-1}}\sum _{i_{j+1}=0}^{N_{j+1}}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}=1,\\ & \mathrm{for}{i}_j=1,\dots N_j\mathrm{and}j=2,\dots k-1,\\ & \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}{z}_{i_{1\cdots }{i}_k}=1,i_k=1,\dots N_k,\\ & z_{i_{1\cdots }{i}_k}\in \left\{0,1\right\},\forall i_n,n=1,\dots ,k,\end{array}& \left(0.4\right)\end{array}$

[0067] where c and z are similarly dimensioned matrices representing costs and hypothetical assignments. Note, in general, for tracking problems these matrices are sparse.

[0068] After formulating the k-dimensional assignment problem as in (0.4), the present invention solves the resulting problem so as to generate the outputs required by devices 110, 112, 122 and 130. For each observation set o_ireceived from converter 104 at time t_iwhere i=1, . . . ,N₁, the computer processes o_iin a batch together with the other observation sets o_i−k+1, . . . ,o_iand the track T_i−k, i.e., T_jis the set of all tracks that have been defined up to but not including o_j. (Note, bold type designations refer to the vector of elements for the indicated time, i.e., the set of all observations in the scan or tracks existing at the time, etc.) The result of this processing is the new set of tracks T_i−k+1and a set of cost weighted possible solutions indicating how the tracks might extend to the current time t_i. At time t_i+1the batch process is repeated using the newest observation set and deleting the oldest. Thus, there is a moving window of observation sets which is shifted forward to always incorporate the most recent observation set. The effect is that input observation sets are reused for (k−1) cycles and then on the observation set's k^threuse each observation within the observation set is integrated into a track.

[0069] FIG. 7 illustrates various processes implemented upon receipt of each observation set. Except for the addition of the k-dimensional assignment solving process 300 and the modification to scoring process 154 to build data structures suitable for process 300, all processes in FIG. 7 are based upon prior art. The following processes 150 and 152 extend previously defined tracks h_i−1based on new observations. Gate formulation and output process 156 determines, for each of the previously defined tracks, a zone wherein the track may potentially extend based on limits of velocity, maneuverability and radar precision. One such technique to accomplish this is the cone method described previously. The definition of the zone is passed to gating process 150. When a new observation set o_iis received, the gating process 150 will match each member observation with the zone for each member of the hypothetical set h_i−1. After all input observations from o_iare processed, the new hypothesis set h_iis generated by extending each track of the prior set of hypothetical tracks h_i−1either with missed detect gap fillers or with all new observation elements satisfying the track's zone. This is a many-to-many matching in that each hypothesis member can be extended to many new observations and each new observation can be used to extend many hypotheses. It, however, is not a full matching in that any hypothesis will neither be matched to all observations nor vice versa. It is this matching characteristic that leads to the sparse matrices involved in the tracking process. Subsequently, gating 150 forwards the new hypothesis set h_ito filtering process 152. Filtering process 152 determines a smooth curve for each member of h_i. Such a smooth curve is more likely than a sharp turn from each point straight to the next point. Further, the filtering process 152 removes small errors that may occur in generating observations. Note that in performing these tasks, the filtering process 152 preferably utilizes a minimization of a least squares test of the points in a track hypothesis or a Kalman filtering approach.

[0070] As noted above, the foregoing track extension process requires knowledge of a previous track. For the initial observations, the following gating process 158 and filtering process 160 determine the “previous track” based on initial observations. In determining the initial tracks, the points from the first observation set form the beginning points of all possible tracks. After observation data from the next observation set is received, sets of simple two-point straight line tracks are defined. Then, promotion, gate formulation, and output step 162 determines a zone in which future extensions are possible. Note that filtering step 160 uses curve fitting techniques to smooth the track extensions depending upon the number of prior observations that have accumulated in each hypothesis. Further note that promotion, gate formulation and output process 162 also determines when sufficient observations have accumulated to form a basis for promoting the track to processes 150 and 152 as described above.

[0071] The output of each of the filtering processes 152 and 160 is a set of hypothetical track extensions. Each such extension contains the hypothetical initial conditions (from the previous track), the list of observations incorporated in the extension, and distance between each observation and the filtered track curve. scoring process 154 determines the figure of merit or cost of an observation being an extension of a track. In one embodiment, the cost is based on the above-mentioned distance although the particular formula for determining the cost is not critical to the present invention. A preferred formula for determining the cost utilizes a negative log likelihood function in which the cost is the negative of the sum of: (a) the logs of the distances normalized by sensor standard deviation parameters, and (b) the log likelihoods for events related to track initiation, track termination, track maneuver, false negatives and false positives. Note that track maneuvers are detected by comparing the previous track curve with the current extension. Further note that some of the other events related to, for example, false negatives and false positives are detected by analyzing the relative relationship of gap fillers in the hypothesis. Thus, after determining that one of these events occurred, a cost for it can be determined based upon suitable statistics tables and system input parameters. The negative log likelihood function is desirable because it permits effective integration of the useful components. Copies of the set of hypothetical track extensions which are scored are subsequently passed directly to one of the gate formulation and output steps 156 and 162. Note that the scoring process 154 also arranges the actual scores in a sparse matrix based upon observation identifiers, and passes them to k-dimensional assignment problem solving process 300.

[0072] The assignment solving process 300 is described below. Its output is simply the list of assignments which constitute the most likely solution of the problem described by Equation (0.4). Note that both gate formulation and output processes 156 and 162 use (at different times) the list of assignments to generate the updated track history T_ito eliminate or prune alternative previous hypotheses that are prohibited by the actual assignments in the list, and, subsequently, to output any required data. Also note that when one of the gate formulation and output processes 156 and 162 accomplish these tasks, the process will subsequently generate and forward the new set of gates for each remaining hypothesis and the processes will then be prepared to receive the next set of observations. In one embodiment, the loop described here will generate zones for a delayed set of observations rather than the subsequent set. This permits processes 156 and 162 to operate on even observation sets while the scoring step 154 and k-dimensional solving process operate on odd sets of observations, or vice versa.

[0073] The assignment solving process 300 permits the present invention to operate with a window size of dimension k−1 for some k≧3. The upper limit on k depends only upon the computational power of the computer 106 and the response time constraints of system 100. The k−1 observation sets within the processing window plus the prior track history result in a k-dimensional Assignment Problem as described by Problem Formulation (0.4). The present invention solves this generalized problem including the processes required to consider false positives and negatives, and also the processes required to consider sparse matrix problem formulations.

[0074] In describing a first embodiment of the k-dimensional assignment solver 300, it is worthwhile to also discuss the process of FIG. 4 which is used by the solver 300. FIG. 4 illustrates use of the Frieze and Yadagar process as shown in prior art for transforming a three-dimensional assignment problem into a two-dimensional assignment problem and then use a hill climbing algorithm to solve the three-dimensional assignment problem. The solver 300 uses a Lagrangian Relaxation technique (well known in the art) to reduce the dimension of an original k-dimensional assignment problem (k>3) down to a three-dimensional problem and then use the process of FIG. 4 to solve the three-dimensional problem. Further note that the Lagrangian Relaxation technique is also utilized by the process of FIG. 4 and that in using this technique the requirement that each point is assigned to one and only one track is relaxed. Instead, an additional cost, which is equal to a respective Lagrangian Coefficient u, is added to the cost or objective function (0.1)(a) whenever a point is assigned to more than one track. This additional cost can be picked to weight the significance of each constraint violation differently, so this additional cost is represented as a vector of coefficients u which are correlated with respective observation points. Hill climbing will then develop a sequence of Lagrangian Coefficients sets designated (u₀, . . . , u_j,u_j+1, . . . ,u_p). That correspond to an optimum solution of the two-dimensional assignment problem. The assignments at this optimum solution are then used to “recover” the assignment solution of the three-dimensional assignment problem.

[0075] In step 200 of FIG. 4, initial values are selected for the u₀coefficients. Because the Lagrangian Relaxation process is iterative, the initial values are not critical and are all initially selected as zero. In step 202, these additional costs are applied to the objective function (0.1)(a). With the addition of the costs u, the goal is still to assign the points which minimize the total cost. This transforms Equation (0.1)(a), written for k=3 and altered to exclude mechanisms related to false positives and negatives, into objective function (1.1)(a). In the first iteration it is not necessary to consider the u matrix because all u values are set to zero. To relax the requirement that each point be assigned to one and only one track, the constraint Equation (0.1)(d) is deleted, thereby permitting points from the last image to be assigned to more than one track. Note that while any axis can be chosen for relaxation, observation constraints are preferably relaxed. The effect of this relaxation is to create a new problem which must have the same solution in the first two axes but which can have a differ solution in the third axis. The result is constraints (1.1)(b-d). $\begin{array}{cc}\begin{array}{ccc}\left(a\right)& \text{Minimize}& \sum _{i_1=0}^{N_1}\sum _{i_2=0}^{N_2}\left\{\left(c_{i_1{i}_2{i}_3}-u_{i_3}\right)\right\}{z}_{i_1{i}_2{i}_3}\\ \left(b\right)& \text{Subject to:}& \sum _{i_2=1}^{N_2}\sum _{i_3=1}^{N_3}{z}_{i_1{i}_2{i}_3}=1,i_1=1,\dots ,N_1,\\ \left(c\right)& & \sum _{i_1=1}^{N_1}\sum _{i_3=1}^{N_3}{z}_{i_1{i}_2{i}_3}\le 1,i_2=1,\dots ,N_2,\\ \left(d\right)& & z_{i_1{i}_2{i}_3}\in \left\{0,1\right\}\forall z_{i_1{i}_2{i}_3}.\end{array}& \left(1.1\right)\end{array}$

[0076] Step 204 then generates from the three-dimensional problem described by Problem Formulation (0.1) a new two-dimensional problem formulation which will have the same solution for the first two indices. As Problem Formulation (1.1) has no constraints on the 3^rdaxis, any value within a particular 3^rdaxis can be used in a solution, but using anything other than the minimum value from any 3^rdaxis has the effect of increasing solution cost. Conceptually, the effect of step 204 is to change the three-dimensional arrays in Problem Formulation (1.1) into two-dimensional arrays as shown in Problem Formulation (1.2) and to generate the new two-dimensional matrix m_iidefined as shown in Equation (1.3). $\begin{array}{cc}\begin{array}{ccc}\left(a\right)& \text{Minimize}& \sum _{i_1=0}^{N_1}\sum _{i_2=0}^{N_2}\left\{\begin{array}{c}\mathrm{Min}:\\ \forall i_3\end{array}\left(c_{i_1{i}_2{i}_3}-u_{i_3}\right)\right\}{z}_{i_1{i}_2}\\ \left(b\right)& \text{Subject to:}& \sum _{i_2=1}^{N_2}\sum _{i_3=1}^{N_3}{z}_{i_2{i}_3}=1,i_1=1,\dots ,N_1,\\ \left(c\right)& & \sum _{i_1=1}^{N_1}\sum _{i_3=1}^{N_3}{z}_{i_1{i}_3}\le 1,i_2=1,\dots ,N_2,\\ \left(d\right)& & z_{i_1{i}_2}\in \left\{0,1\right\}\forall z_{i_1{i}_2},\end{array}& \left(1.2\right)\end{array}$  m_ii=Minimum{c_ii−u_j|i=1, . . . ,N_k}.  (1.3)

[0077] The cost or objective function for the reduced problem as defined by (1.2)(a), if evaluated at all possible values of u is a surface over the domain of u. This surface is referred to as Φ(u) and is non-smooth but provable convex (i.e., it has a single peak and several other critical characteristics which form terraces). Due to the convex characteristics of Φ(u), the results from solving Problem Formulation (1.2) at any particular u_jcan be used to generate a new set of coefficients u_j+1whose corresponding cost value is closer to the peak of Φ(u). The particular set of Lagrangian Coefficients that will generate the two-dimensional problem resulting in the maximum cost is designated u_p. To recover the solution to the three-dimensional assignment problem requires solving the Equation (1.2) problem corresponding to u_p.

[0078] In step 206, the two-dimensional problem is solved directly using a technique known to those skilled in the art such as Reverse Auction for the corresponding cost and solution values. This is the evaluation of one point on the surface or for the first iteration Φ(u₀).

[0079] Thus, after this first iteration, the points have been assigned based on all “u” values being arbitrarily set to zero. Because the “u” values have been arbitrarily assigned, it is unlikely that these assignments are correct and it is likely that further iterations are required to properly assign the points. Step 208 determines whether the points have been properly assigned after the first iteration by determining if for this set of assignments whether a different set of “u” values could result in a higher total cost. Thus, step 208 is implemented by determining the gradient of objective function (1.2)(a) with respect to u_j. If the gradient is substantially non-zero (greater than a predetermined limit) then the assignments are not at or near enough to the peak of the surface Φ(u)(decision 210), and the new set of Lagrangian Coefficients u_j+1is determined.

[0080] Hill climbing Step 212 determines the u_j+1values based upon the u_jvalues, the direction resulting from protecting the previous gradient into the U³domain, and a step size. The solution value of the two-dimensional problem is the set of coefficients that minimize the two-dimensional problem and the actual cost at the minimum. Those coefficients augmented by the coefficients stored in m_iipermit the evaluation (but not the minimization) of the cost term in Problem Formulation (1.1). These two cost terms are lower and upper bounds on the actual minimized cost of the three-dimensional problem, and the difference between them in combination with the gradient is used to compute the step size.

[0081] With this new set of u values, steps 202-210 are repeated as a second iteration. Steps 212 and 202-210 are repeated until the gradient as a function of u determined in step 208 is less than the predetermined limit. This indicates that the u_pvalues which locate the peak area of the Φ(u) surface are determined and that the corresponding Problem Formulation (1.2) has been solved. Step 214 will attempt to use the assignments that resulted from this particular two-dimensional assignment problem to recover the solution of the three-dimensional assignment problem as described below. If the limit was chosen properly so that the u values are close enough to the peak, this recovery will yield the proper set of assignments that rigidly satisfies the constraint that each point be assigned to one and only one track. However, if the u values are not close enough to the peak, then the limit value for decision 210 is reduced and the repetition of steps 212 and 202-210 is continued.

[0082] Step 214 recovers the three-dimensional assignment solution by using the assignment values determined on the last iteration through step 208. Consider the two-dimensional z assignment matrix to have 1's in the locations specified by the list L₁=(a_i, b_i)_i=1^N. If the three-dimensional z matrix is specified by placing 1's at the location indicated by the list L₂=(a_i, b_i, m_ab)_i=1^Nthen the result is a solution of Problem Formulation (1.1). Let L₃=(m_ab)_i=1^Nbe the list formed by the third index. If each member of L₃is unique then the L₂solution satisfies the third constraint so it is a solution to Problem Formulation (0.1). When this is not the case, recovery determines the minimal substitutions required within list L₃so that it plus L₁will be a feasible solution, i.e., a solution which satisfies the constraints of a problem formulation, but which may not optimize the objective function of the problem formulation. This stage of the recovery process is formulated as a two-dimensional assignment problem as follows. Form a new cost matrix [c_i,j]_i,j=1^Nwhere c_i,j=c_ab,jfor j=1 . . . N_iand the N_iterm is the total number of cost elements in the selected row of the three-dimensional cost matrix. Attempt to solve this two-dimensional problem for the minimum using two constraints sets. If a feasible solution is found then the result will have the same form as list L₁. Replace the first set of indices by the indicated (a_i, b_i) pairs taken from list L₁and the result will be a feasible solution of Problem Formulation (0.4). If no feasible solution to the new two-dimensional problem exists then further effort to locate the peak of Φ(u) is required.

[0083] I.1. Generalization to a Multi-Dimensional Assignment Solving Process

[0084] Let M be a fixed integer and assume that M is the dimension of the initial assignment problem to be solved. Thus, initially, the result of the scoring step 154 is a M-dimensional Cost Matrix which is structured as a sparse matrix (i.e., only a small percentage of the entries in the cost and assignment matrices are filled or non-zero). Individual cost elements represent the likelihood that a track T_ias extended by the set of observations {O_ij|i=1, . . . ,M−1}, is not valid. Because the matrix is sparse the list of cost elements is stored as a packed list, and then for each dimension of the matrix, a vector of a variable length list of pointers to the cost elements is generated and stored. This organization means that for a particular observation o_ijfor the j^thlist in the i^thvector will be a list of pointers to all hypotheses in which O_ijparticipates. This structure is further explained in the following section dealing with problem partitioning.

[0085] The objective of the assignment solving process is to select from the set of all possible combinations of track extensions a subset that satisfies two criteria. First, each point in the subset of combinations should be assigned to one and only one track and therefore, included in one and only one combination of the subset, and second, the total of the scoring sums for the combinations of the subset should be minimized. This yields the following M-dimensional equations where k=M: $\begin{array}{cc}\begin{array}{c}\begin{array}{ccc}\left(a\right)& \text{Minimize}{v}_k\left(z^k\right)=\sum _{i_1=0}^{N_1}\cdots \sum _{i_k=0}^{N_k}{c}_{i_{1\cdots }{i}_k}^k{z}_{i_{1\cdots }{i}_k}^k& \end{array}\\ \begin{array}{ccc}\left(b\right)& \text{Subject to:}& \sum _{i_2=0}^{N_2}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^k=1,\\ \left(c\right)& & \sum _{i_1=0}^{N_1}\cdots \sum _{i_{j-1}=0}^{N_{j-1}}\sum _{i_{j+1}=0}^{N_{j+1}}\dots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^k=1,\\ & & \mathrm{for}{i}_j=1,\dots ,N_j\mathrm{and}j=2,\dots ,k-1,\\ \left(d\right)& & \sum _{i_1=0}^{N_1}\cdots \sum _{i_{k-1}=0}^{N_{k-1}}{z}_{i_{1\cdots }{i}_k}^k=1,i_k=1,\dots ,N_k,\\ \left(e\right)& & z_{i_{1\cdots }{i}_k}^k\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1,\dots ,l_k,\end{array}\end{array}& \left(1.4\right)\end{array}$

[0086] and where c^kis the cost matrix [c_{i. . . i}^k] which is a function of the distance between the observed point z^kand the smoothed track determined by the filtering step, and v_kis the cost function. This set of equations is similar to the set presented in Problem Formulation (0.4) except that it includes the subscript and superscript k notation. Thus, in solving the M-dimensional Assignment Problem the invention reduces this problem to an (M−1)-dimensional Assignment Problem and then to an (N−2)-dimensional Assignment Problem, etc. Further, the symbol kε{3, . . . ,M} customizes Problem Formulation (1.4) to a particular relaxation level. That is, the notation is used to reference data from levels relatively removed as in c^k+lare the cost coefficients which existed prior to this level of relaxed coefficients c_k. Note that actual observations are numbered from 1 to N_i, where N_iis the number of observations in observation set i. Further note that the added O observation in each set of observations is the unconstrained “gap filler.” This element serves as a filler in substituting for missed detects, and a sequence of these elements including only one true observation represents the possibility that the observation is a false positive. Also note that by being unconstrained a gap filler may be used in as many hypotheses as required.

[0087] While direct solution to (1.4) would give the precise assignment, the solution of k-dimensional equations directly for large k is too complex and time consuming for practice. Thus, the present invention solves this problem indirectly.

[0088] The following is a short description of many aspects of the present invention and includes some steps according to the prior art. The first step in solving the problem indirectly is to reduce the complexity of the problem by the previously known and discussed Lagrangian Relaxation technique. According to the Lagrangian Relaxation technique, the absolute requirement that each point is assigned to one and only one track is relaxed such that for some one image, points can be assigned to more than one track. However, a penalty based on a respective Lagrangian Coefficient u^kis added to the cost function when a point in the image is assigned to more than one track. The Lagrangian Relaxation technique reduces the complexity or “dimension” of the formulation of the assignment problem because constraints on one observation set are relaxed. Thus, the Lagrangian Relaxation is performed iteratively to repeatedly reduce the dimension until a two-dimensional penalized cost function problem results as in Problem Formulation (1.1). This two-dimensional problem is solved then directly by a previously known technique such as Reverse Auction. The penalized cost function for the two-dimensional problem defines a valley or convex shaped surface which is a function of various sets of {u^k|k=3, . . . ,M} penalty values and one set of assignments for the points in two dimensions. That is, for each particular u³there is a corresponding two-dimensional penalized cost function problem and its solution. Note that the solution of the two-dimensional penalized cost function problem identifies the set of assignments for the particular u³values that minimize the penalized cost function. However, these assignments are not likely to be optimum for any higher dimensional problem because they were based on an initial arbitrary set of u_kvalues. Therefore, the next step is to determine the optimum assignments for the related three-dimensional penalized cost function problem. There exists a two-dimensional hill shaped function Φ which is a graph of the minimums of all penalized cost functions at various sets of assignments in two dimensions. For the three-dimensional problem, the function Φ can be defined based on the solution to the foregoing two-dimensional penalized cost function. By using the current u³values and the {u^k|k>3} values originally assigned, the gradient of the hill-shaped function Φ is determined, which points toward the peak of the hill. By using the gradient and u³values previously selected for the one point on the hill (corresponding to the minimum of the penalized cost function Φ) as a starting point, the u³values can be found for which the corresponding problem will result in the peak of the function Φ. The solution of the corresponding two-dimensional problem is the proper values for two of the three sets of indices in the three-dimensional problem. These solution indices can select a subsection of the cost array which maps to a two-dimensional array. The set of indices which minimize the two-dimensional assignment problem based on that array corresponds to the proper assignment of points in the third dimension. The foregoing “recovery” process was known in the prior art, but it is modified here to adjust for the sparse matrix characteristic. The next task is to recover the solution of the proper u⁴values for the four-dimensional problem. The foregoing hill climbing process will not work again because the foregoing hill climbing process when required to locate the four-dimensional u⁴values for the peak of Φ³requires the exact definition of the function Φ³(as was available in the case of Φ²) or an always less than approximation of Φ³, whereas the iteration can result in a greater than approximation of Φ³. According to the present invention, another three-dimensional function Ψ is defined which is a “less than approximation” the three-dimensional function Φ and which can be defined based on the solution to the two-dimensional penalized cost function and the previously assigned and determined u^kvalues. Next, the gradient of the function Ψ is found and hill climbing technique used to determine the u⁴values at the peak. Each selected u⁴set results in a new three-dimensional problem and requires the two-dimensional hill climbing based upon new two-dimensional problems. At the peak of the three-dimensional function Ψ, the solution values are a subset of the values required for the four-dimensional solution. Recovery processing extends the subset to a complete solution. This process is repeated iteratively until the u^kvalues that result in a corresponding solution at the peak of the highest order functions Ψ and Φ are found. The final recovery process then results in the solution of k-dimensional problem.

[0089] FIG. 5 illustrates process 300 in more detail.

[0090] I.2. Problem Formulation

[0091] In problem formulation step 310, all data structures for the subsequent steps are allocated and linked by pointers as required for execution efficiency. The incoming problem is partitioned as described in the subsequent section. This partitioning has the effect of dividing the incoming problem into a set of independent problems and thus reducing the total workload. The partitioning process depends only on the actual cost matrix so the partitioning can and is performed for all levels of the relaxation process.

[0092] I.2.1. Relaxation and Recovery

[0093] Step 320 begins the Lagrangian Relaxation process for reducing the M-dimensional problem by selecting all Lagrangian Coefficient u^Mpenalty values initially equal to zero. The Lagrangian Coefficients associated with the M^thconstraint set are an (N_M+1)-element vector. The reduction of this M-dimensional problem in step 322 to a (M−1)-dimensional problem uses the two step process described above. First, a penalty based on the value of the respective u^Mcoefficient is added to the cost function when a point is assigned to more than one track and then the resultant cost function is minimized. However, during the first iteration, the penalty is zero because all u^Mvalues are initially set to zero. Second, the requirement that no point from any image can be assigned to more than one track is relaxed for one of the images. In the extreme this would allow a point from the relaxed image to be associate with every track. However, the effect of the previous penalty would probably mean that such an association would not minimize the cost. The effect of the two steps in combination is to remove a hard constraint while adding the penalty to the cost function so that it operates like a soft constraint. For step 322 this two-step process results in the following penalized cost function problem with k=M: $\begin{array}{cc}\begin{array}{c}\begin{array}{cc}\left(a\right)& {\Phi }_k\left(u^k\right)\equiv \text{Minimize:}{\phi }_k\left(u^k,z^k\right)\\ & =\sum _{i_1=0}^{N_1}\dots \sum _{i_k=0}^{N_k}{c}_{i_{1\cdots }{i}_k}^k{z}_{i_{1\cdots }{i}_k}^k-\\ & \sum _{i_k=0}^{N_k}{u}_{i_k}^k\lfloor \sum _{i_1=0}^{N_1}\dots {\sum _i}_{k-1}^{N_{k-1}}{z}_{i_{1\cdots }{i}_k}^k-1\rfloor \\ & =\sum _{i_1=0}^{N_1}\dots \sum _{i_k=0}^{N_k}\left(c_{i_{1\cdots }{i}_k}^k-u_{i_k}^k\right)z_{i_{1\cdots }{i}_k}^k+\sum _{i_k=0}^{N_k}{u}_{i_k}^k\end{array}\\ \begin{array}{ccc}\left(b\right)& \text{Subject to:}& \sum _{i_2=0}^{N_2}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^k=1,i_1=1,\dots N_1,\\ \left(c\right)& & \sum _{i_1=0}^{N_1}\cdots \sum _{i_{j-1}=0}^{N_{j-1}}\sum _{i_{j+1}=0}^{N_{j+1}}\dots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^k=1\\ & & \mathrm{for}{i}_j=1,\dots N_j\\ & & \mathrm{and}j=2,\dots k-1,\\ \left(d\right)& & z_{i_{1\cdots }{i}_k}^k\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1\dots l_k.\end{array}\end{array}& \left(1.5\right)\end{array}$

[0094] Because the constraint on single assignment of elements from the last image has been eliminated, an (M−1)-dimensional problem can be developed by eliminating some of the possible assignments. As shown in Equations (1.6), this is done by selecting the smallest cost element from each of the M^thaxis vectors of the cost matrix. Reduction in this manner yields a new, lower-order penalized cost function defined by Equations (1.6) which has the same minimum cost as does the objective function defined by (1.5)(a) above.

[0095] The cost vectors are selected as follows. Define a new cost array c_{i. . . i}^k−1by: $\begin{array}{cc}\begin{array}{c}{c}_{i_1{\mathrm{\dots i}}_k}^k=\mathrm{Minimum}\left\{c_{i_1{\mathrm{\dots i}}_k}^k-u_{i_k}^k❘i_k=0,1,\dots ,N_k\right\},\\ \mathrm{for}\left(i_1,\dots ,i_{k-1}\right)\ne \left(0,\dots ,0\right),\\ c_{0\mathrm{\dots 0}}^{k-1}=\sum _{i_k=0}^{N_k}\min \left\{0,c_{0\mathrm{\dots 0}}^k-u_{i_k}^k\right\}.\end{array}& \left(1.6\right)\end{array}$

[0096] The resulting (M−1)-dimensional problem is (where k=M): $\begin{array}{cc}\begin{array}{c}{\Phi }_k\left(u^k\right)=\text{Minimize:}{v}_{k-1}\left(z^{k-1}\right)\\ =\sum _{i_1=0}^{N_1}\cdots \sum _{i_{k-1}=0}^{N_{k-1}}{C}_{i_{1\cdots }{i}_{k-1}}^{k-1}{z}_{i_{1\cdots }{i}_{k-1}}^{k-1}\end{array}& \left(1.7\right)\\ \begin{array}{cc}\text{Subject to:}& \sum _{i_2=0}^{N_2}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^{k-1}=1,i_1=1,\dots N_1,\\ & \sum _{i_1=0}^{N_1}\cdots \sum _{i_{j-1}=0}^{N_{j-1}}\sum _{i_{j+1}=0}^{N_{j+1}}\cdots \sum _{i_k=0}^{N_k}{z}_{i_{1\cdots }{i}_k}^{k-1}=1,\\ & \mathrm{for}{i}_j=1,\dots N_j\mathrm{and}j=2,\dots k-1,\\ & \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-2}=0}^{N_{k-2}}{z}_{i_{1\cdots }{i}_{k-1}}^{k-1}=1,i_{k-1}=1,\dots N_{k-1},\\ & z_{i_{1\cdots }{i}_k}^k\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1\dots l_{k-1}.\end{array}& \end{array}$

[0097] Assignment Problem (1.4) and Problem Formulation (1.7) differ only in the dimension M vs. M−1, respectively. An optimal solution to (1.7) is also an optimal solution to equation (1.5). This relationship is the basis for an iterative sequence of reductions indicated by steps 320-332 through 330-332 and 200-204 in which the penalized cost function problem is reduced to a two-dimensional problem. As these formula will be used to describe the processing at all levels, the lowercase k is used except where specific reference to the top level is needed. In step 206, the two-dimensional penalized cost function is solved directly by the prior art Reverse Auction technique. Each execution of 206 produces two results, the set of z²values that minimize the problem and the cost that results v²when these z values are substituted into the objective function (1.7)(a).

[0098] In step 208, according to the prior art, solution z²values are substituted into the two-dimensional derivative of the surface Φ₂. The result indicates how the value of u³should be adjusted so as to perform the hill climbing function. As was previously described the objective is to produce a sequence of u_i³values which ends when the u_p³value is in the domain of the peak of the surface Φ. The section “Determining Effective Gradient” describes how new values are computed and how it is determined that the current u_i^kpoints to the peak of Φ₂. When no further adjustment is required the flow moves to step 214 which will attempt to recover the three-dimensional solution as previously described. When further adjustment is required then the flow progresses to step 212 and the new values of u^kare computed. At the two-dimensional level the method of the prior art could be used for the hill climbing procedure. However, it is not practical to use this prior art hill climbing technique to determine the updated Lagrangian Coefficients u^kor the Max on the next (or any) higher order surface Φ because the next (or any) higher dimensional function Φ cannot be defined based on known information.

[0099] Instead, the present invention defines a new function based on known information which is useful for hill climbing from the third to the fourth and above dimensions, i.e., u^kvalues which result in z values that are closer to the proper z values for the highest k-dimension. This hill climbing process (which is different than that used in the prior art of FIG. 4 for recovering only the three-dimensional solution) is used iteratively at all lower dimensions k of the M-dimensional problem (including the three-dimensional level where it replaces prior art) even when k is much larger than three. FIG. 6 helps to explain this new hill climbing technique and illustrates the original k-dimensional cost function v_kof Problem Formulation (1.4). However, the actual k-dimensional cost surface v_kmdefined by (1.4)(a) comprises scalar values at each point described by k-dimensional vectors and as such can not be drawn. Nevertheless, for illustration purposes only, FIG. 6 ignores the reality of ordering vectors and illustrates a concave function v_k(z^k) to represent Equation (1.4). The surface v_kis illustrated as being smooth to simplify the explanation although actually it can be imagined to be terraced. The goal of the assignment problem is to find the values of {overscore (z)}^k; these values minimize the k-dimensional cost function v_k.

[0100] For purposes of explanation, assume that in FIG. 6, k=4 (the procedure is used for all k≧3). This problem is reduced by two iterations of Lagrangian Relaxation to a two-dimensional penalized cost function φ_j^(k−1). This cost function, and all other cost functions described below, are also non-smooth and continuous but are illustrated in FIG. 6 as smooth for explanation purposes. Solving the φ_j^(k−1)problem results in one set of z²assignments and the value of Φ_(k−1)at the point u_ij². A series of functions φ_j^(k−1), . . . ,φ₁^(k−1)each generated from a different u³are shown. The particular series illustrates the process of locating the peak of Φ_(k−1). The two-dimensional penalized cost functions φ_j^(k−1), . . . ,φ₁^(k−1)can be solved directly. Each such solution provides the information required to calculate the next u³value. Each iteration of the hill climbing improves the selection of u³values, i.e., yields φ²problem whose solution is closer to those at the solution of the φ³problem. The result of solving φ₁^(k−1)is values that are on both Φ_(k−1)and Φ_k. FIG. 6 illustrates the surface Φ_kwhich comprises the minimums of all k-dimensional penalized cost function surfaces φ_k, i.e., if the Problem Formulations (1.5) and (1.7) were solved at all possible values of u^kthe function Φ_k(u_k) would result. The surface Φ_kis always less than the surface v_kexcept at the peak as described in Equation (1.8) and its maximum occurs where the surface v_kis minimum. Because the surface Φ_krepresents the minimum of the surface Φ_k, any point on the surface Φ_kcan mapped to the surface v_k. The function Φ_kprovides a lower bound on the minimization problem described by (1.4)(a). Let {overscore (z)}^kbe the unknown solution to Problem Formulation (1.4) and note that:

[0101] Φ^k(u^k)≦v_k({overscore (z)}^k)≦v_k−1(z^k).  (1.8)

[0102] Consequently, the values z^kat the peak of Φ_k(i.e., the greatest lower bound on the cost of the relaxed problem), can be substituted into the k-dimensional penalized cost function to determine the proper assignments. Consequently, the present invention attempts to find the maximum on the surface Φ_k. However, it is not possible to use the prior art hill climbing to hill climb to the peak of Φ_kbecause the definition of Φ_krequires exact knowledge of lower order functions Φ. As the solution of φ_l^kis not the exact solution, in that higher order values of u are not yet optimized, its value can be larger than the true peak of Φ_k. As such it is not a lower bound on v_kand it can not be used to recover the required solution.

[0103] Instead, the present invention defines all auxiliary functions Ψ_kwhich are based only on the solution to the two-dimensional penalized cost function problem, lower order values z^kand values u^kdetermined previously by the reduction process. The function Ψ_kis a less than approximation of Φ_k, and its gradient is used for hill climbing to its peak. The values z^kat the peak of the function Ψ_kare then substituted into Problem Formulation (1.5) to determine the proper assignments. To define the function Ψ_k, the present invention explicitly makes the function Φ_k(u^k) a function of all higher order sets of Lagrangian Coefficients with the expanded notation: Φ_k(u^k;u^k+1, . . . ,u^k) . Then, a new set of functions Ψ, is defined recursively, using the Φ_k's domain: $\begin{array}{cc}{\psi }_3\left(u^3\right)=v_2+\sum _{i_3=0}^{N_3}{u}_{i_3}^3={\Phi }_3\left(u^3;u^4,\dots ,u^K\right),& \left(1.9\right)\end{array}$

[0104] where v₂is the solution value for the most recent two-dimensional penalized cost function problem. For k>3 $\begin{array}{cc}{\psi }_k\left(u^3,\dots ,u^{k-1};u^k\right)=\{\begin{array}{cc}{\Phi }_k\left(u^k;u^{k+1},\dots ,u^K\right),& \mathrm{if}\mathrm{known},\\ {\psi }_{k-1}\left(u^3,\dots ,u^{k-2};u^{k-1}\right)+\sum _{i_k=0}^{N_k}{u}_{i_k}^k,& \mathrm{otherwise}.\end{array}& \left(1.10\right)\end{array}$

[0105] From the definitions of Φ_kand v_k(Problem Formulation (1.7) compared with Problem Formulation (1.5)): ${\Phi }_k\left(u^k;u^{k+1},\dots ,u^M\right)=v_{k-1}\left(z^{k-1}\right)+\sum _{i_k=0}^{N_k}{u}_{i_k}^k,$

[0106] it follows that: ${\psi }_3\left(u^3\right)=v_2+\sum _{i_3=0}^{N_3}{u}_{i_3}^3={\Phi }_3\left(u^3;u^4,\dots ,u^M\right)\le v_3^3\left(z^3\right)$

[0107] and with that Equation (1.8) is extended to:

[0108] Ψ_k(u³, . . . ,u^k−1;u^k)≦Φ_k(u^k;u^k+1, . . . ,u^M)≦v_k({overscore (u)}^k)≦v_k(u^k)  (1.11)

[0109] This relationship means that either Φ_kor Ψ_kmay be used in hill climbing to update the Lagrangian Coefficients u. Φ_kis the preferred choice, however it is only available when the solution to Problem Formulation (1.5) is a feasible solution to Problem Formulation (1.4) (as in hill climbing from the second to third dimension which is why prior art worked). For simplicity in implementation, the function Ψ_kis defined so that it equals Φ_kwhen either function could be used. It is therefore always used, even for hill climbing from the second dimension. The use of the function Ψ which is based on previously assigned or determined u values and not higher order u values which are not yet determined, is an important feature of the present invention.

[0110] I.2.2. Determining Effective Gradient

[0111] After the function Ψ_kis defined, the next steps of hill climbing/peak detection are to determine the gradient of the function Ψ_k, determine an increasing portion of the gradient and then move up the surface Ψ_kin the direction of this increasing portion of the gradient. As shown in FIG. 6, any upward step on the surface Ψ_k, for example, to the minimum of the φ^k_jwill yield a new set of values u^k(to the “left”) that is closer to the ideal set of values u^kwhich correspond to the minimum of the function v_k. While it is possible to iteratively step up this surface Ψ_kwith steps of fixed size and then determine if the peak has been reached, the present invention optimizes this process by determining the single step size from the starting point at the minimum of φ^k_ithat will jump to the peak and then calculate the values u^kat the peak. Once the values u at the peak of Ψ_kare determined, then the values u^kcan be substituted into Problem Formulation (1.5) to determine the proper assignment. (However, in a more realistic example, where k is much greater than three, then the values u^kat the peak of the function Ψ along with the lower-order values u^kand those assigned and yielded by the reduction steps are used to define the next higher level function Ψ. This definition of a higher order function Ψ and hill climbing process are repeated iteratively until Ψ_k, the peak of Ψ_k, and the values u^kat the peak of Ψ_kare identified.) The following is a description of how to determine the gradient of each surface Ψ and how to determine the single step size to jump to the peak from the starting point on each surface Ψ.

[0112] As noted above, each surface Ψ is non-smooth. Therefore, if a gradient was taken at a single point, the gradient may not point toward the peak. Therefore, several gradients (a “bundle”) are determined at several points on the surface Ψ_kin the region of the starting point (i.e., minimum of φ_k) and then averaged. Statistically, the averaged result should point toward the peak of the surface Ψ_k. Wolfe's Conjugate Subgradient Algorithm (P. Wolfe. A method of conjugate subgradients for minimizing non-differentiable functions. Mathematical Programming Study,3:145-173, 1975; P. Wolfe. Finding the nearest point in a polytope. Mathematical Programming Study,11;128-149, 1976) for minimization was previously known in another environment to determine a gradient of a non-smooth surface using multiple subgradients and can be used with modification in the present invention. The modification to Wolfe's algorithm uses the information generated for Ψ_k(u^k3, . . . ,u^kk−2;u^kk−1) as the basis for calculating the subgradients. The definition of a subgradient v of Ψ_k(u^k3, . . . ,u^kk) is any member of the subdifferential set defined as:

[0113] ∂Ψ_k(u)={vεR^N+1|(Ψ_k(u^k3, . . . ,u^kk−1;u′)−Ψ_k(u^k3, . . . ,u^k))≧v^T(u′−u^k)∀u′εR^N+1},

[0114] where v^Tis the transpose of v.

[0115] Next, a subgradient vector is determined from this function. If z^kis the solution of Problem Formulation (1.5), then differentiating Ψ_k(u³₁, . . . ,u^k−1;u^k) with respect to u_ik^kand evaluating the result with respect to the current selection matrix z^kyields a subgradient vector: $g=\left(0,\left(1-\sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_k}{z}_{i_1{\mathrm{\dots i}}_k}^k|i_k=1,\dots ,N_k\right)\right).$

[0116] The iterative nature of the solution process at each dimension yields a set of such subgradients. Except for a situation described below, where the resultant averaged gradient does not point toward the peak, the most recent set of such subgradients are saved and used as the “bundle” for the peak finding process for this dimension. For example, at the level k there is a bundle of subgradients of the surface Ψ_knear the minimum of the surface φ^k_idetermined as a result of solving Problem Formulation (1.4) at all lower levels. This bundle can be averaged to approximate the gradient. Alternately, the previous bundle can be discarded so as to use the new value to initiate a new bundle. This choice provides a way to adjust the process to differing classes of problems, i.e., when data is being derived from two sensors and the relaxation proceeds from data derived from one sensor to the other then the prior relaxation data for the first sensor could be detrimental to performance on the second sensors data.

[0117] I.2.3. Step Size and Termination Criterion

[0118] After the average gradient of the surface Ψ_kis determined, the next step is to determine a single step that will jump to the peak of the surface Ψ_k. The basic strategy is first to specify an arbitrary step size, and then calculate the value of Ψ_kat this step size in the direction of the gradient. If the value of Ψ_kis larger than the previous one, this probably means that the step has not yet reached the peak. Consequently, the step size is doubled and a new value of Ψ_kis determined and compared to the previous one. This process is repeated until the new value of Ψ_kis less than the previous one. At that time, the step has gone too far and crossed over the peak. Consequently, the last doubling is rolled-back and that step size is used for the iteration. If the initial estimated Ψ_kis less than previous value then the step size is decreased by 50%. If this still results in a smaller value of Ψ_k, then the last step is rolled back and the previous step size is decreased by 25%. The following is a more detailed description of this process.

[0119] With a suitable bundle of subgradients determined as just described, Wolfe's algorithm can be used to determine the effective subgradient d and the upgraded value u_j+1^k. From the previous iteration, or from an initial condition, there exists a step length value t. The value,

[0120] u₊=u_j^k+td

[0121] is calculated as an estimate of u_j+1^k. To determine if the current step size is valid we evaluate Ψ_k(u³, . . . ,u^k−2;u₊) . If the result represents an improvement then we double the step size. Otherwise we halve the step size. In either case a new u₊ is calculated. The doubling or halving continues until the step becomes too large to improve the result, or until it becomes small enough to not degrade the result. The resulting suitable step size is saved with d as part of the subgradient bundle. The last acceptable u₊ is assigned to u_j+1^k.

[0122] Three distinct criteria are used to determine when u_j^kis close enough to {overscore (u)}^k:

[0123] 1. The Wolfe's algorithm criterion of d=0 given that the test has been repeated with the bundle containing only the most recent subgradient.

[0124] 2. The difference between the lower bound Φ_k(u^k) and the upper bound v_k(z^k,u^k) being less than a preset relative threshold. (Six percent was found to be an effective threshold for radar tracking problems.)

[0125] 3. An iteration count being exceeded.

[0126] The use of limits on the iteration are particularly effective for iterations at the level 3 through (n−1) as these iterations will be repeated so as to resolve higher order coefficient sets. With limited iterations the process is in general robust enough to improve the estimate of upper order Lagrangian Coefficients. By limiting the iteration counts then the total processing time for the algorithm becomes deterministic. That characteristic means the process can be effective for real time problems such as radar, where the results of the last scan of the environment must be processed prior to the next scan being received.

[0127] I.2.4. Solution Recovery

[0128] The following process determines if the u^kvalues at what is believed to be the peak of the function Ψ_kadequately approximate the u^kvalues at the peak of the corresponding function Φ_k. This is done by determining if the corresponding penalized cost function is minimized. Thus, the u^kvalues at the peak of the function Ψ_kare first substituted into Problem Formulation (1.5) to determine a set of z assignments for k−1 points. During the foregoing reduction process, Problem Formulation (1.7) yielded a tentative list of k selected z points that can be described by their indices as: {(i₁^j. . . i_k−1^j)}_j=1^N, where N₀is the number of cost elements selected into the solution. One possible solution of Problem Formulation (1.4) is the solution of Problem Formulation (1.5) which is described as {(i₁^j. . . i_k−1^jm _i^k. . . i_k−1)}_j=1^Nwith m_i^k. . . i_k−1as it 5 was defined in Problem Formulation (1.6). If this solution satisfies the k^thconstraint set, then it is the optimal solution.

[0129] However, if the solution for Problem Formulation (1.5) is not feasible (decision 355), then the following adjustment process determines if a solution exists which satisfies the k^thconstraint while retaining the assignments made in solving Problem Formulation (1.7). To do this, a two-dimensional cost matrix is defined based upon all observations from the k^thset which could be used to extend the relaxed solution. $\begin{array}{cc}{h}_{\mathrm{jl}}=c_{i_1,\dots ,i_{k-1},l}^k& \left(1.12\right)\\ \mathrm{for}l=0,\dots ,N_k& \\ \mathrm{and}j=0,\dots ,N_0.& \end{array}$

[0130] If the resulting two-dimensional assignment problem, $\begin{array}{cc}\mathrm{Minimize}\sum _{j=0}^{N_0}\sum _{l=0}^{N_1}{h}_{\mathrm{jl}}{w}_{\mathrm{jl}}& \left(1.13\right)\\ \text{Subject to:}& \\ \begin{array}{cc}\sum _{l=0}^{N_k}{w}_{\mathrm{jl}}=1& j=1,\dots ,N_0,\\ \sum _{l=0}^{N_k}{w}_{\mathrm{jl}}=1& j=1,\dots ,N_0,\\ w_{\mathrm{jl}}\in \left\{0,1\right\}& j=0,\dots ,N_0,l=0,\dots ,N_k\end{array}& \end{array}$

[0131] has a feasible solution, then the indices of that solution map to the solution of Problem Formulation (1.4) for the k-dimensional problem. The first index in each resultant solution entry is the pointer back to an element of the {(i₁^j. . . i_k−1^jm_i. . . i_k−1)}_j=1^Nlist. That element supplies the first k−1 indices of the solution. The second index of the solution to the recovery problem is the k^thindex of the solution. Together these indexes specify the values of z^kthat solve Problem Formulation (1.4) at the k^thlevel.

[0132] If Problem Formulation (1.13) does not have a feasible solution, then the value of u_i^kwhich was thought to represent u_p^kis not representative of the actual peak and further iteration at the k^thlevel is required. This decision represent the branch path from step 214 and equivalent steps.

[0133] I.3. Partitioning

[0134] A partitioning process is used to divide the cost matrix that results from the Scoring step 154 into as many independent problems as possible prior to beginning the relaxation solution. The partitioning process is included with the Problem Formulation Step 310. The result of partitioning is a set of problems to be solved, i.e., there will be p₁problems that consist of a single hypothesis, p₂problems that consist of two hypothesis, etc. Each such problem is a group in that one or more observations or tracks are shared between members of the group.

[0135] In partitioning to groups no consideration is given to the actual cost values. The analysis depends strictly on the basis of two or more cost elements sharing the same specific axis of the cost matrix. In a two-dimensional case two cost elements must be in the same group if they share a row or a column. If the two elements are in the same row, then each other cost element that is also in the same row, as well as any cost elements that are in columns occupied by members of the row must be included in the group. The process continues recursively. In literature it is referred to as “Constructing a Spanning Forest.” The k-dimensional case is analogous to the two-dimensional case. The specific method we have incorporated is a depth first search, presented by Aho, Hopecraft, and Ullman (A. V. Aho, J. E. Hopcroft, and J. D. Ullman. Design and Analysis of Computer Algorithms. Addison-Wesley, MA, 1974).

[0136] The result of partitioning at level M is the set of problems described as {P_ij|i=1, . . . ,p_jand j=1, . . . ,N}, where N is the total number of hypothesis. The problems labeled {P_i1|i=1, . . .p₁} are the cases where there is only one choice for the next observation at each scan and that observation could be used for no other track, i.e., it is a single isolated track. The hypothesis must be included in the solution set and no further processing is required.

[0137] As hypothesis are constructed the first element is used to refer to the track ID. Any problem in the partitioned set which does not have shared costs in the first scan represent a case where a track could be extended in several ways but none of the ways share an observation with any other track. The required solution hypothesis for this case is the particular hypothesis with the maximum likelihood. For this case all likelihoods are determined as was described in Scoring and the maximum is selected.

[0138] In addition to partitioning at the level M, partitioning is applied to each subsequent level M−1 . . . ,2. For each problem that was not eliminated by the prior selection, partitioning is repeated, ignoring observations that are shared in the last set of observations. Partitioning recognizes that relaxation will eliminate the last constraint set and thus partitioning is feasible for the lower level problems that will result from relaxation. This process is repeated for all levels down to k=3. The full set of partitionings can be performed in the Problem Formulation Step 310, prior to initiating the actual relaxation steps. The actual two-dimensional solver used in step 206 includes an equivalent process so no advantage would be gained by partitioning at the k=2 level.

[0139] There are two possible solution methods for the remaining problems. “Branch and Bound” as was previously described, or the relaxation method that this invention describes. If any of the partitions have 5-10 possible tracks and less than 50 to 20 hypotheses, then the prior art “Branch and Bound” algorithm generally executes faster than does the relaxation due to its reduced level of startup overhead. The “Branch and Bound” algorithm is executed against all remaining M level problems that satisfy the size constraint. For the remainder the Relaxation algorithm is used. The scheduling done in Problem Formulation allows each Problem Formulation (1.5) cost matrix resulting from the first step of relaxation to be partitioned. The resulting partitions can be solved by any of the four methods: isolated track direct inclusion, isolated track tree evaluation, small group “Branch and Bound” or an additional stage of relaxation as has been fully described.

[0140] The partitioning after each level of Lagrangian Relaxation is effective because when a problem is relaxed, the constraint that requires each point to be assigned to only one track is eliminated (for one image at a time) . Therefore, two tracks previously linked by contending for one point will be unlinked by the relaxation which permits the point to be assigned to both tracks. The fruitfulness of partitioning increases for lower and lower levels of relaxation.

[0141] The following is a more detailed description of the partitioning method. Its application at all but the level M depends upon the relaxation process described in this invention. The recursive partitioning is therefore a distinct part of this invention. The advantage of this method is greatly enhanced by the sparse cost matrix resulting from tracking problems. However the sparse nature of the problem requires special storage and search techniques.

[0142] A hypothesis is the combination of the cost element c_nthe selection variable z_nand all observations that made up the potential track extension. It can be written as, h_n={c_n, z_n, {o_n=o_k|k=1, . . . ,M_iε{1 . . . ,N_k}}}, i.e., cost, selection variable and observations in the hypothesis. Here nε{1, . . .N}, where N is the total number of hypotheses in the problem. While the cost and assignment matrices were previously referenced, these matrices are not effective storage mechanisms for tracking applications. Instead the list of all hypothesis and sets of lists for each dimension that reference the hypothesis set are stored. The hypothetical set in list form is:

[0143] {i h_n}_n=1^n=N.

[0144] For each dimension k=1, . . . ,M there exists a set of lists, with each list element being a pointer to a particular hypothesis:

[0145] $L_{k_i}={\left\{p_{k_j}\text{❘}{k}_j=1,\dots ,N_{k_i}\right\}}_{i=1,k=1}^{i=N_{k_i},k=M},$

[0146] where N_kis a number of hypothesis containing the i^thobservation from scan k. This structure is a multiply linked list in that any observation is associated with a set of pointer to all hypothesis it participates in, and any hypothesis has a set of pointers to all observations that formed it. (These pointers can be implemented as true pointers or indices depending upon the particular programming language utilized.)

[0147] Given this description of the matrix storage technique then the partitioning technique is as follows: Mark the first list L_kand follow out all pointers in that list to the indicated hypothesis $h_{{\mathrm{Pk}}_i}\mathrm{for}i=1,\dots ,N_{k_i}.$

[0148] Mark all located hypothesis, and for each follow pointers back the particular L_kfor k=1 , . . . ,M. Those L's if not previously marked get marked and also followed out to hypothesis elements and back to L's. When all such L's or h's being located are marked, then an isolated sub-problem has been identified. All marked elements can be removed from the lists and stored as a unique problem to be solved. The partitioning problem then continues by again starting at the first residual set L. When none remain, the original problem is completely partitioned.

[0149] Isolated problems can result from one source track having multiple possible extensions or from a set of source tracks contending for some observations. Because one of the indices of k (in our implementation it is k=1) indicates the source track then it is possible to categorize the two problem types by observing if the isolated problem includes more than one L-list from the index level associated with tracks.

[0150] I.4. Subsequent Track Prediction

[0151] As noted above, after the points are assigned to the respective tracks, some action is taken such as controlling take-offs and landings to avoid collision, advising airborne aircraft to change course, warning airborne aircraft of an adjacent aircraft in a commercial environment, or aiming and firing an anti-aircraft (or anti-missile) missile, rocket or projectile, or taking evasive action in a military environment. Also, the tracking can be used to position a robot to work on an object. For some or all of these applications, usually the tracks which have just been identified are extended or extrapolated to predict a subsequent position of the aircraft or other object. The extrapolation can be done in a variety of prior art ways; for example, straight line extensions of the most likely track extension hypothesis, parametric quadratic extension of the x versus time, y versus time, etc. functions that are the result of the filtering process described earlier, least square path fitting or Kalman filtering of just the selected hypothesis.

[0152] The process of gating, filtering, and gate generation as they were previously described require that a curve be fit through observations such that fit of observations to the likely path can be scored and further so that the hypothetical target future location can be calculated for the next stage of gating. In the implementation existing, quadratic functions have been fit through the measurements that occur within the window. A set of quadratics is used, one per sensor measurement. To calculate intercept locations these quadratics can be converted directly into path functions. Intercept times are then calculated by prior methods based upon simultaneous solution of path equations.

[0153] The use of the fitted quadratics is not as precise as more conventional filters like the Kalman Filter. They are however much faster and sufficiently accurate for the relative scores required for the Assignment problem. When better location prediction is required then the assignment problem is executed to select the solution hypothesis and based upon the observations in those hypothesis the more extensive Kalman filter is executed. The result is tremendous computation savings when compared with the Kalman Filter being run on all hypothetical tracks.

[0154] Based on the foregoing, apparatus and methods have been disclosed for tracking objects. However, numerous modifications and substitutions can be made without deviating from the scope of the present invention. For example, the foregoing functions Ψ can also be defined as recursive approximation problem in which several values of higher order u^kvalues are used to eliminate the higher than approximation characteristic of the function Φ_k. The hill climbing of the function Ψ can be implemented by a high order hill climbing using the enhanced function Φ_k. Although the result would not be as efficient it seems likely that the method would converge. Therefore, the invention has been disclosed by way of example and not limitation, and reference should be made to the following claims to determine the scope of the present invention.

[0155] The following description provides an alternative second embodiment of the multi-dimensional assignment solving process of section I.1. However, in order to clearly describe the present alternative embodiment, further discussion is first presented regarding the data assignment problem of partitioning measurements into tracks and false alarms and, in particular, the representation of multi-dimensional assignment problems.

[0156] II.1. Formulation of the Assignment Problem

[0157] The goal of this section is to explain the formulation of the data association problems, and more particularly multi-dimensional assignment problems, that govern large classes of problems in centralized or hybrid centralized-sensor level multisensor/multitarget tracking. The presentation is brief; technical details are presented for both track initiation and maintenance in (A. B. Poore. Multidimensional assignment formulation of data association problems arising from multitarget tracking and multisensor data fusion. Computational Optimization and Applications,3:27-57, 1994) for non-maneuvering targets and in (A. B. Poore. Multidimensional assignments and multitarget tracking: Partitioning data sets. In P. Hansen, I. J. Cox, and B. Julesz, editors, DIMACS Series in Discrete Mathematics and Theoretical Computer Science, volume 19, pages 169-198, Providence, R.I., 1995. American Mathematical Society) for maneuvering targets. Note that the present formulation can also be modified to include target features (e.g., size and type) into the scoring process 154.

[0158] The data assignment problems for multisensor and multitarget tracking considered here are generally posed as that of maximizing the posterior probability of the surveillance region (given the data) according to $\begin{array}{cc}\mathrm{Maximize}\left\{\frac{P\left(\Gamma =\gamma |Z^M\right)}{P\left(\Gamma ={\gamma }^0|Z^M\right)}|\gamma \in {\Gamma }^{*}\right\},& \left(2.1\right)\end{array}$

[0159] where Z^Mrepresents M data sets, γ is a partition of indices of the data (and thus induces a partition of the data), Γ* is the finite collection of all such partitions, Γ is a discrete random element defined on Γ*, γ⁰is a reference partition, and P(Γ=γ|Z^M) is the posterior probability of a partition γ being true given the data Z^M. The term partition is defined below.

[0160] Consider M observation sets Z(k) (k=1, . . . , N) each of N_kobservations {z_i^k}_i=1^Nand let Z^Mdenote the cumulative data set defined by

[0161] Z(k)={z_i^k}_i=1^Nand Z^m={Z(1), . . . , Z(M)},  (2.22)

[0162] respectively. In multisensor data fusion and multitarget tracking the data sets Z(k) may represent different classes of objects, and each data set can arise from different sensors. For track initiation the objects are measurements that must be partitioned into tracks and false alarms. In the formulation of track extensions a moving window over time of observations sets is used. The observation sets will be measurements which are: (a) assigned to existing tracks, (b) designated as false measurements, or (c) used for initiating tracks. However, note that alternative data objects instead of observation sets may also be fused such as in sensor level tracking wherein each sensor forms tracks from its own measurements and then the tracks from the sensors are fused in a central location. Note that, as one skilled in the art will appreciate, both embodiments of the present invention may be used for this type of data fusion.

[0163] The data assignment problem considered presently is represented as a case of set partitioning defined in the following way. First, for notational convenience in representing tracks, a zero index is added to each of the index sets in (2.2), a gap filler z₀^kis added to each of the data sets z(k) in (2.2), and a “track of data” is defined as (z_i¹, . . . , z_i^N), where i_kcan now assume the values of 0. A partition of the data will refer to a collection of tracks of data or track extensions wherein each observation occurs exactly once in one of the tracks of data and such that all data is used up. Note that the occurrence of the gap filler is unrestricted. The gap filler z₀^kserves several purposes in the representation of missing data, false observations, initiating tracks, and terminating tracks. Note that a “reference partition” may be defined which is a partition in which all observations are declared to be false.

[0164] Next under appropriate independence assumptions it can be shown that: $\begin{array}{cc}\frac{P\left(\Gamma =\gamma |Z^M\right)}{P\left(\Gamma ={\gamma }^0|Z^M\right)}\equiv L_{\gamma }\equiv \prod _{\left(i_1{\mathrm{\dots i}}_N\right)\in \gamma }{L}_{i_1{\mathrm{\dots i}}_M},& \left(2.3\right)\end{array}$

[0165] where L_{i. . . i}^kis a likelihood ratio containing probabilities for detection, maneuvers, and termination as well as probability density functions for measurement errors, track initiation and termination. Then if c_{i. . . i}^k=−ln(L_{i. . . i}^k), $\begin{array}{cc}-\ln \lfloor \frac{P\left(\Gamma =\gamma |Z^M\right)}{P\left(\Gamma ={\gamma }^0|Z^M\right)}\rfloor =\sum _{\left(i_1,\dots ,i_M\right)\in \gamma }{c}_{i_1{\mathrm{\dots i}}_N}.& \left(2.4\right)\end{array}$

[0166] Using (2.3) and the zero-one variable z_{i. . . i}^k=1, if (i_l, . . . ,i_M) ε γ, and 0, otherwise, one can then write the problem (2.1) as the following M-dimensional assignment problem (as also presented in section Problem Formulation (1.4) with k=M): $\begin{array}{cc}\begin{array}{ccc}\left(a\right)& \mathrm{Minimize}& \sum _{i_1=0}^{N_1}\dots \sum _{i_M=0}^{N_M}{C}_{i_1\dots i_M}{Z}_{i_1\dots i_M}\\ \left(b\right)& \text{Subject to:}& \sum _{i_2=0}^{N_2}\dots \sum _{i_M=0}^{N_M}{z}_{i_1\dots i_M}=1,i_1=1,\dots ,N_1,\\ \left(c\right)& & \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}\sum _{i_{k+1}=0}^{N_{k+1}}\dots \sum _{i_M=0}^{N_M}{z}_{i_1\dots i_M}=1,\\ & & \mathrm{for}{i}_k=1,\dots ,N_k\mathrm{and}k=2,\dots ,M-1,\\ \left(d\right)& & \sum _{i_1=0}^{N_1}\dots \sum _{i_{M-1}=0}^{N_{M-1}}{z}_{i_1\dots i_M}=1,i_M=1,\dots ,N_M,\\ \left(e\right)& & z_{i_1\dots i_M}\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1,\dots ,i_M,\end{array}& \left(2.5\right)\end{array}$

[0167] where c₀. . . ₀is arbitrarily defined to be zero. Here, each group of sums in th e constraints (2.5)(b)-(2.5)(e) represents the fact that each non-gap filler observation occurs exactly once in a “track of data.” One can modify this formulation to include multi-assignments of one, some, or all the actual observations. The assignment problem (2.5) is changed accordingly. For example, if z_i^kis to be assigned no more than, exactly, or no less than n_i^ktimes, then the “=1” is changed to ≦, =, ≧n_i^k, respectively.

[0168] Expressions for the likelihood ratios L_{i. . . i}^kappearing in the costs c_{i. . . i}=−ln (L_{i. . . i}) are well-known. In particular, discussions of these ratios may be found in (A.B. Poore. Multidimensional assignment formulation of data association problems arising from multitarget tracking and multisensor data fusion, Computational Optimization and Applications,3:27-57, 1994; A.B. Poore. Multidimensional assignments and multitarget tracking: Partitioning data sets. In P. Hansen, I. J. Cox, B. Julesz, editors, DIMACS Series in Discrete Mathematics and Theoretical Computer Science, volume 19, pages 169-198, Providence, R.I., 1995. American Mathematical Society). Furthermore, the likelihood ratios are easily modified to include target features and to account for different sensor types. Also note that in track initiation, the M observation sets provide observations from M sensors, possibly all the same. Additionally note that for track maintenance, a sliding window of M observation sets and one data set containing established tracks may be used. In this latter case, the formulation is the same as above except that the dimension of the assignment problem is now M+1.

[0169] II.2. Overview of the New Lagrangian Relaxation Algorithms

[0170] Having formulated an M-dimensional assignment problem (2.5), we now turn to a description of the Lagrangian relaxation algorithm. Subsection II.2.1 below presents many of the relaxation properties associated with the relaxation of an n-dimensional assignment problem to an m-dimensional one via a Lagrangian relaxation of n−m sets of constraints wherein m<n≦M and preferably in the present invention embodiment n−m>1. Although any n−m constraint sets can be relaxed, the description here is based on relaxing the last n−m sets of constraints and keeping the first m sets. Given either an optimal or suboptimal solution of the relaxed problem, a technique for recovering a feasible solution of the n-dimensional problem is presented in subsection II.2.2 below and an overview of the Lagrangian relaxation algorithm is given in subsection II.2.3.

[0171] The following notation will be used throughout the remainder of the work. Let M be an integer such that M≧3 and let nε{3, . . . , M}. The n-dimensional assignment problem is $\begin{array}{cc}\begin{array}{ccc}\left(a\right)& \text{Minimize:}{v}_n\left(z\right)=& \sum _{i_1=0}^{N_1}\dots \sum _{i_n=0}^{N_n}{c}_{i_1\dots i_n}^n{z}_{i_1\dots i_n}^n\\ \left(b\right)& \text{Subject to}& \sum _{i_2=0}^{N_2}\dots \sum _{i_n=0}^{N_n}{z}_{i_1\dots i_n}^n=1,i_1=1,\dots ,N_1,\\ \left(c\right)& & \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}\sum _{i_{k+1}=0}^{N_{k+1}}\dots \sum _{i_n=0}^{N_n}{z}_{i_1\dots i_n}^n=1,\\ & & \mathrm{for}{i}_k=1,\dots ,N_k\mathrm{and}k=2,\dots ,n-1,\\ \left(d\right)& & \sum _{i_1=0}^{N_1}\dots \sum _{i_{n-1}=0}^{N_{n-1}}{z}_{i_1\dots i_n}^n=1,i_n=1,\dots ,N_n,\\ \left(e\right)& & z_{i_1{\mathrm{\dots i}}_n}^n\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1,\dots ,i_n.\end{array}& \left(2.6\right)\end{array}$

[0172] To ensure that a feasible solution of (2.6) always exists, all variables with exactly one nonzero index (i.e., variables of the form z₀ⁿ. . . _0i0. . . ₀ⁿfor i_k≠0) are assumed free to be assigned and the corresponding cost coefficients are well-defined.

[0173] II.2.1. The Lagrangian Relaxed Assignment Problem

[0174] The n-dimensional assignment problem (2.6) has n sets of constraints. A (N_k+1)-dimensional multiplier vector associated with the k^thconstraint set will be denoted by u^k=(u₀^k, u₁^k, . . . , u_M^k)^Twith u₀^k=0 and k=1 , . . . , n. The n-dimensional assignment problem (2.6) is relaxed to an m-dimensional assignment problem by incorporating n−m of the n sets of constraints into the objective function (2.6)(a). Although any constraint set can be relaxed, the description of the relaxation procedure for (2.6) will be based on the relaxation of the last n−m sets of constraints. The relaxed problem is $\begin{array}{cc}\begin{array}{c}{\phi }_{\mathrm{mn}}\left(u^{m+1},\dots ,u^n\right)\equiv \mathrm{Minimize}{\phi }_{\mathrm{mn}}\left(z^n;u^{m+1},\dots ,u^n\right)\\ \equiv \sum _{i_1=0}^{N_1}\dots \sum _{i_n=0}^{N_n}{c}_{i_1\dots i_n}^n{z}_{i_1\dots i_n}^n+\sum _{k=m+1}^n\sum _{i_k=0}^{N_k}\\ u_{i_k}^k\left[\sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}\sum _{i_{k+1}=0}^{N_{k+1}}\dots \sum _{i_n=0}^{N_n}{z}_{i_1\dots i_n}^n-1\right]\\ \equiv \sum _{i_1=0}^{N_1}\dots \sum _{i_n=0}^{N_n}\left[c_{i_1\dots i_n}^n+\sum _{k=m+1}^n{u}_{i_k}^k\right]z_{i_1\dots i_n}^n-\\ \sum _{k=m+1}^n\sum _{i_k=0}^{N_k}{u}_{i_k}^k\\ \text{Subject To}\sum _{i_2=0}^{N_2}\dots \sum _{i_n=0}^{N_n}{z}_{i_1\dots i_n}^n=1,i_1=1,\dots ,N_1,\\ \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}\sum _{i_{k+1}=0}^{N_{k+1}}\dots \sum _{i_n=0}^{N_n}{z}_{i_1\dots i_n}^n=1\\ \mathrm{for}i=1,\dots ,N_k\mathrm{and}k=2,\dots ,m,\\ z_{i_1\dots i_n}^n\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1,\dots ,i_n,\end{array}& \left(2.7\right)\end{array}$

[0175] wherein we have added the multiplier u₀^k≡0 for notational convenience. Thus, the multiplier u^kε^N+1with u₀^k=0 is fixed.

[0176] An optimal (or suboptimal) solution of (2.7) can be constructed from that of an m-dimensional assignment problem. To show this, define for each (i_l, . . . , i_m) an index (j_m+1, . . . ,j_n)=(j_m+1(i_l, . . . , i_m), . . . , j_n(i₁, . . . , i_m)) and a new cost function c_{i. . . i}^mby $\begin{array}{cc}\left(j_{m+1},\dots ,j_n\right)=\arg \min \left\{c_{i_1\dots i_n}^n+\sum _{k=m+1}^n{u}_{i_k}^k|i_k=0,\dots ,N_k\mathrm{and}k=m+1,\dots ,n\right\},& \left(2.8\right)\\ c_{i_1\dots i_m}^m=c_{i_1\dots i_m{j}_{m+1}\dots j_n}+\sum _{k=m+1}^n{u}_{j_k}^k\mathrm{for}\left(i_1,\dots ,i_m\right)\ne \left(0,\dots 0\right),c_{0\dots 0}^m=\sum _{i_{m+1}}^{N_{m+1}}\dots \sum _{i_n=0}^{N_0}\mathrm{Minimum}\left\{0,c_{0\dots 0i_{m+1}\dots i_n}^n+\sum _{k=m+1}^n{u}_{j_k}^k\right\}.& \end{array}$

[0177] If (j_m+1, . . . , j_n) is not unique, choose the smallest such j_m+1, amongst those (n−m)-tuples with the same j_m+1choose the smallest j_m+2, etc., so that (j_m+1, . . . , j_n) is uniquely defined.) Using the cost coefficients defined in this way, the following m-dimensional assignment problem is obtained: $\begin{array}{cc}\begin{array}{c}{\hat{\phi }}_{\mathrm{mn}}\left(u^{m+1},\dots ,u^n\right)=\mathrm{Minimize}{\hat{\phi }}_{\mathrm{mn}}\left(z^m;u^{m+1},\dots ,u^n\right)\\ \equiv v_m\left(z^m\right)\equiv \sum _{i_1=0}^{N_1}\dots \sum _{i_m=0}^{N_0}{c}_{i_1\dots i_m}^m{z}_{i_1\dots i_m}^m\end{array}& \left(2.9\right)\\ \begin{array}{c}\text{Subject to:}\sum _{i_2=0}^{N_1}\dots \sum _{i_m=0}^{N_m}{z}_{i_1\dots i_m}^m=1,i_1=1,\dots ,N_1,\\ \sum _{i_1=0}^{N_1}\dots \sum _{i_{k-1}=0}^{N_{k-1}}\sum _{i_{k+1}=0}^{N_{k+1}}\dots \sum _{i_m=0}^{N_0}{z}_{i_1\dots i_m}^m=1,\\ \mathrm{for}{i}_k=1,\dots ,N_k\mathrm{and}k=2,\dots ,m-1,\\ \sum _{i_1=0}^{N_1}\dots \sum _{i_{m-1}=0}^{N_{m-1}}{z}_{i_1\dots i_m}^m=1,i_m=1,\dots ,N_m,\\ z_{i_1\dots i_m}^m\in \left\{0,1\right\}\mathrm{for}\mathrm{all}{i}_1,\dots ,i_m.\end{array}& \end{array}$