SYSTEM AND METHOD FOR MODELING ADVANCED AUTOMOTIVE SAFETY SYSTEMS

Abstract

A system and methods are disclosed for providing an integrated software development environment for the design, verification, and validation of advanced automotive safety systems. The system allows automotive software to be developed on a host computer using a collection of computer programs running simultaneously as processes and synchronized by a central process. The software disclosed uses separate synchronized processes, permitting signals from disparate sources to be generated by a simulation running on the host computer or from actual sensors and data bus signals coming from and going to actual vehicle hardware which is connected to their bus counterparts in the host computer on a real-time basis. The methods provide a data model that first extends the capabilities of the physical data model and then translates, gates, optimizes, fuses, filters, and manages the physical representation of the logical model into a state estimation of the situation around the vehicle.

Description

This application further incorporates by reference U.S. patents:

U.S. Pat. No. 6,629,033 Issued Sep. 30, 2003 Titled—OPEN COMMUNICATION SYSTEM FOR REAL-TIME MULTIPROCESSOR APPLICATIONS;

U.S. Pat. No. 6,771,208 Issued Aug. 3, 2004 Titled—MULTI SENSOR SYSTEM;

U.S. Pat. No. 7,146,260 Issued Dec. 5, 2006 Titled—METHOD AND APPARATUS FOR DYNAMIC CONFIGURATION OF MULTIPROCESSOR SYSTEM;

U.S. Pat. No. 6,792,351 Issued Sep. 14, 2004 Titled—METHOD AND APPARATUS FOR MULTI-VEHICLE COMMUNICATION;

U.S. Pat. No. 6,778,073 Issued Aug. 17, 2004 Titled—METHOD AND APPARATUS FOR MANAGING AUDIO DEVICES;

U.S. Pat. No. 6,615,137 Issued Sep. 2, 2003 Titled—METHOD AND APPARATUS FOR TRANSFERRING INFORMATION BETWEEN VEHICLES;

U.S. Pat. No. 7,178,049 Issued Feb. 13, 2007 Titled—METHOD FOR MULTI-TASKING MULTIPLE JAVA VIRTUAL MACHINES;

U.S. Pat. No. 7,337,650 Issued Mar. 4, 2008 Titled—SYSTEM AND METHOD FOR ALIGNING SENSORS ON A VEHICLE.

COPYRIGHT NOTICE

Contained herein is material that is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure, as it appears in the United States Patent and Trademark Office patent file or records, but otherwise reserves all rights to the copyright whatsoever. The following notice applies to the software, screenshots and data as described below and in the drawings hereto and All Rights Reserved.

TECHNICAL FIELD OF THE INVENTION

This invention relates generally to statistical modeling and analysis of automotive safety systems, sensors and application requirements for both performance and cost tradeoffs. This invention relates specifically to the process of designing, building, testing and analyzing system data in a hierarchy of statistical models, driven by the need to implement advanced safety features for the transportation industry.

BACKGROUND

Next generation automotive systems will include numerous electronic systems such as; passive safety systems for airbag deployment and anti-skid braking features; active safety for collision warning and collision avoidance features; and convenience features for Blind Spot Detection (BSD) and Adaptive Cruise Control (ACC). Today the transportation industry is rapidly moving toward solutions that support and enable these systems and new features with the ultimate goal of profitably and developing a more survivable vehicle at costs their customers are willing to pay.

These new systems will require information about the environment, targets in the environment and our relationship to them. This information will be generated from multiple sources to include data from a new class of sensor, termed non-contact sensor, and include non-contact sensors like radar, optical, laser, ultrasonic, etc. These non-contact sensors will generate data about the environment including range measurements and target classification. These systems are similar to those found on advanced tactical fighter aircraft today.

The process of manufacturing a vehicle is not much different today then is was 50 years ago; features are specified for a desired model and year; parts designed; competitively purchased; assembled; and shipped as a complete vehicle. If another feature is desired, it is integrated into the design and bolted on as a complete feature; this approach is referred to as “federated” systems. A key tenant to low cost systems in the future will be a deliberate movement away from federated features and systems, to fully integrated systems design and integrated systems.

Now with a focus on active safety systems, the sensors required to support the features identified above are to some extent common; an optical sensor that supports one optical application could be used for a second optical application; two sensors could be used for one application. In either case, simply sharing data from sensors or other sources to improve knowledge of the environment is incorrectly termed “sensor-fusion” or “multi sensor-fusion” by the transportation industry suppliers.

Those skilled in the art of state estimation, robotics, advanced defense avionics understand academically that sensor-fusion is the art of combining sensory data or data derived from disparate sources such that the resulting information is in some sense “better” than would be possible when these sources were used individually. This process is predicated on the covariance (or the measure of how much two variables vary together) of non-independent sources. The term “better” in the case above can mean more accurate, more complete, more dependable, or refer to the result of an emerging view or state estimation.

The data sources for a fusion process are not specified to originate from identical sources or sensors which may or may not be spatially and temporally aligned. Further one can distinguish direct fusion, indirect fusion and fusion of the outputs of the former two. Direct fusion is the fusion of sensor data from a set of heterogeneous or homogeneous sensors, soft sensors, and history values of sensor data, while indirect fusion uses information sources like a prior knowledge about the environment and human input. Sensor fusion is also known as “multi-sensor data fusion” and is a subset of information fusion through an implementation of the probability theory.

Probability theory is the mathematical study of phenomena characterized by randomness or uncertainty. More precisely, probability is used for modeling situations when the result of a measurement, realized under the same circumstances, produces different results.

Mathematicians and actuaries think of probabilities as numbers in the closed interval from 0 to 1 assigned to “events” whose occurrence or failure to occur is random. Two crucial concepts in the theory of probability are those of a random variable and of the probability distribution of a random variable.

Implementing the features described above with affordable instruments requires reliable real-time estimates of system state. Unfortunately, the complete state is not always observable. State Estimation takes all the data obtained and uses it to determine the underlying behavior of the system at any point in time. It includes fault detection, isolation and continuous system state estimation.

There are two parts to state estimation: modeling and algorithms. The overall approach is to use a model to predict the behavior of the system in a particular state, and then compare that behavior with the actual measurements from the instruments to determine which state or states is the most likely to produce the observed system behavior.

This is not well understood or currently implemented in the transportation industry; the approach understood and practiced is logical decisions in linear and deterministic systems. If use cases require higher confidences in measurements, instrument specifications are tightened resulting in the undesired effect of cost and schedule increases. The environment we live and operate in is neither linear nor deterministic; use cases are infinite; and the perverse variability of the targets and potential maneuvers cannot be modeled. The variability of the problem identified above includes aspects other than just spatial (i.e. range and bearing to a target); temporal relationships are part of the fundamental intellectual structure (together with space and number) within which events must be sequenced, quantify the duration of events, quantify the intervals between them, and compare the kinematics of objects.

The automotive industry today implements features in such a way that all aspects are contained within the system (federated) and therefore reasonably controlled. Sharing information like target reports with other features is anticipated and desired; however measurements reported with respect to the integrated system will be historical in nature and received asynchronously. Timing errors can induce more error in the system than a bad measurement. These and other issues can be addressed with the introduction of a suite of modeling tools based on re-thinking the approach of federated systems and focusing on an integrated systems solution based on state estimation.

Central to the successful implementation of the advanced safety systems discussed above is ensuring ability for the system to cope with and recover from emergency situations. If one or more emergency conditions arise, systems of the future must quickly initiate and successfully execute procedures to mitigate the condition and then recover; these procedures must be performed under tight timing constraints, e.g., pre-crash systems.

This patent describes a system and methods necessary to implement a design methodology that will facilitate and support advanced safety systems design, test, verification, and validation; with emphasis on reliable fault tolerant operation. The implementation of this system and methods is termed “Algorithm Prototyping Analysis and Test” or APAT.

As it stands today, there is much art published documenting the research and development in the area of procedure analysis and design. However, there are no systematic and rigorous methodologies for designing procedures to be used in advanced safety systems for the transportation industry. This is currently viewed as a serious shortcoming, since these high-risk and complex systems employ procedures and action sequences that can and do impact life or death results.

In the future it will be the responsibility of the onboard computers to automatically initiate and execute procedures and recovery sequences. Therefore modeling, analysis, verification, and design of these highly critical methods, procedures and recovery sequences are required and are thus the focus of this patent.

In a discussion of Prior Art, U.S. Pat. No. 7,006,950, Statistical Modeling and Performance Characterization of a Real-time Dual Camera Surveillance System, Greiffenhagen, et al.; the invention relates to a method for visually detecting and tracking an object through a space. The method derives statistical models for errors, including quantifying an indexing step performed by an indexing module, and tuning system parameters. Further the method applies a likelihood model for candidate hypothesis evaluation and object parameters estimation for locating the object. This invention relates specifically to computer vision systems, more particularly to a system having computationally efficient real-time object detection, tracking, and zooming capabilities. The need arises for methods of modeling more than 1 data source against features that require an infinitely variable combination of instruments and measurements.

In a discussion of Prior Art, U.S. Pat. No. 6,028,548, Vehicle Collision Radar with Randomized FSK Waveform, Farmer; describes an automotive radar for improved immunity to jamming from other automotive radars utilizing common modeling methods such as Auto Regressive Modeling (ARMA) and Minimum Variance Spectral Estimation which are just two such methods that would be applicable. It is recognized in the art that modeling sequences implemented as tools or embedded systems will yield desired results. As with U.S. Pat. No. 7,006,950, the modeling sequences identified are based on the improvement of federated devices.

In a discussion of Prior Art, U.S. Pat. No. 7,079,993, Automated Generator of Optimal Models for the Statistical Analysis of Data, Stephenson, et al; describes an automated process for producing accurate statistical models from sample data tables designed to solve for optimal parameters of each statistical model considered, computes test statistics and degrees of freedom in the model, and uses these test statistics and degrees of freedom to establish a complete ordering of the statistical models. Stephenson further describes how the process arrives at a statistical model that is highest in the ordering and is thus deemed most suitable to model the sample data table. These techniques are basically described in the general art of data or information fusion and modeling. It would not be obvious that general and well published statistical techniques should be applied to the art of automotive safety systems design.

Therefore, a need exists for a system and method for modeling sensor system inputs in a synthetic environment against truth tables to achieve optimal performance and cost parameters for the design and implementation of advanced automotive safety systems.

SUMMARY OF THE INVENTION

Although the best understanding of the present invention will be had from a through reading of the specification and claims presented below, this summary is provided in order to acquaint the reader with some of the new and useful features of the present invention. Of course, this summary is not intended to be a complete litany of all of the features of the present invention, nor is it intended in any way to limit the breadth of the claims, which are presented at the end of the detailed description of this application.

In the preferred embodiment, a system and methods are disclosed for providing an Algorithm Prototyping, Analysis and Test (APAT) through an integrated framework for dynamic data modeling and application development. The system provides a data model that first extends the capabilities of the physical data model and then translates, gates, optimizes, fuses, filters and manages the physical representation of the logical model into a state estimation of the situation around the vehicle; this is known as Situation Awareness (SA); the implementation of this core is called SACore (Situation Awareness Core). And because the application framework is integrated with SACore, once the user configures the data model, the features and benefits of the application framework become automatically available to the developer as part of the decision process for cost-performance optimization.

The data-modeling framework of APAT enables a user to define a data model or version of SACore by using a series of graphical user interfaces (GUI). APAT and SACore provide functionality that physical databases cannot provide by dynamically translating the asynchronous data generated into a corresponding state estimation. Users can utilize the GUI to modify the data model by adjusting expected error, rate or core operation, and the system output is automatically adjusted to accommodate the changes. The APAT system and methods shields both the application developer and the sensor supplier from the complexities, limitations and costs of single sensor systems.

Therefore methods need to be developed to allow automotive system designers the ability to characterize numerous sensor types and sensor suppliers within a specific type as to the actual capability of the sensor in terms of measured variances in the reports. Second the methods need to allow the designers the freedom to place the sensor on the vehicle with other sensors of like type or different types; these may have coincident, mostly coincident, or non-coincident fields of view. Further these methods need to allow the designer the ability to specify tolerances of the dynamic system, such as target maneuver expectations. It would be desirable to the system designer if methods could be presented that allowed all system variables to be extracted into a single file allowing the designer the flexibility of changing any variable in a complex system of variables and structured filters without the tedious need of hand coding the changes. In an n-scalable system, wherein the number of sensors, and number of sensor reports are variable, the resulting identity matrix or unit matrix of size n is represented by the n-by-n square matrix, where n-by-n matrices are used to represent linear transformations from an n-dimensional vector space to itself. Defining and coding a 2-by 2 square matrix then changing the square matrix to a 3-by-3 or n-by-n is tedious and probably beyond the skill set of the automotive system designer, therefore methods are presented as the objects of the invention that allow rapid prototyping and auto-code generation for the development, modeling, implementation and testing of advanced automotive safety systems.

It is an objective of the invention to describe methods for rapidly defining or redefining an automotive multi-sensor data fusion system comprising at least one of a vehicle, sensors, a system model of expected behaviors of targets with consideration for system variables, system error and measurement tolerances. The underlying system is based on an n-scalable, asynchronous, measurement oriented, measurement to track, variable structure, interacting multiple model estimation filter, utilizing a nearest neighbor joint probabilistic data association filter using optimal routines, with two dynamic switchable data association filters in place of the optimal JPDA using sub-optimal routines including one of cheap JPDA and sub-optimal JPDA; the system is implemented with consideration for out-of-sequence measurements, a track management system and a reporting system based on polar coordinates and a method for executing a debiased conversion, translation and rotation into a uniform Cartesian Coordinate system.

Another object of the invention is to define the modeling framework including identifying the system variables by group, wherein the first group comprises sensor parameters and includes such parameters as the number of sensors in the system; the relationship of the sensor frames to the vehicle reference; the expected report limits; report outputs; 1-sigma covariance terms.

Another object of the invention is to define the second group of modeling variables as the target maneuver accelerations that will be acceptable in measurement reports, wherein the modeling environment is at least two dimensional, and wherein at least three models are anticipated and they comprise near constant velocity, cross range acceleration, and down range acceleration, and the values are in terms of expected acceleration g's of maneuver and these accelerations are introduced into the filter structure as noise values.

Another object of the invention is to describe a method for reducing processing complexity, when desired, by setting the number of filters implemented on any pairings of sensors, by identifying number of possibilities, with subsequent values set as noise for each option chosen.

Another object of the invention is to describe the filter setup comprising at least one of housekeeping in nature; values set that define the function of the model values that set the starting uncertainty covariance terms expected from each sensor with respect to the dimensional states identified which further defined down range and cross range values for position, acceleration and velocity.

Other features of the present invention will be apparent from the accompanying drawings and from the detailed description that follows.

Aspects and applications of the invention presented here are described below in the drawings and detailed description of the invention. Unless specifically noted, it is intended that the words and phrases in the specification and the claims be given their plain, ordinary, and accustomed meaning to those of ordinary skill in the applicable arts. The inventors are fully aware that they can be their own lexicographers if desired. The inventors expressly elect, as their own lexicographers, to use only the plain and ordinary meaning of terms in the specification and claims unless they clearly state otherwise and then further, expressly set forth the “special” definition of that term and explain how it differs from the plain and ordinary meaning. Absent such clear statements of intent to apply a “special” definition, it is the inventors' intent and desire that the simple, plain and ordinary meaning to the terms be applied to the interpretation of the specification and claims.

The inventors are also aware of the normal precepts of English grammar. Thus, if a noun, term, or phrase is intended to be further characterized, specified, or narrowed in some way, then such noun, term, or phrase will expressly include additional adjectives, descriptive terms, or other modifiers in accordance with the normal precepts of English grammar. Absent the use of such adjectives, descriptive terms, or modifiers, it is the intent that such nouns, terms, or phrases be given their plain, and ordinary English meaning to those skilled in the applicable arts as set forth above.

Further, the inventors are fully informed of the standards and application of the special provisions of 35 U.S.C. §112, 6. Thus, the use of the words “function,” “means” or “step” in the Detailed Description or Description of the Drawings or claims is not intended to somehow indicate a desire to invoke the special provisions of 35 U.S.C. §112, 6, to define the invention. To the contrary, if the provisions of 35 U.S.C. §112, 6 are sought to be invoked to define the inventions, the claims will specifically and expressly state the exact phrases “means for” or “step for, and will also recite the word “function” (i.e., will state “means for performing the function of [insert function]”), without also reciting in such phrases any structure, material or act in support of the function. Thus, even when the claims recite a “means for performing the function of . . . ” or “step for performing the function of . . . ”, if the claims also recite any structure, material or acts in support of that means or step, or that perform the recited function, then it is the clear intention of the inventors not to invoke the provisions of 35 U.S.C. §112, 6. Moreover, even if the provisions of 35 U.S.C. §112, 6 are invoked to define the claimed inventions, it is intended that the inventions not be limited only to the specific structure, material or acts that are described in the preferred embodiments, but in addition, include any and all structures, materials or acts that perform the claimed function as described in alternative embodiments or forms of the invention, or that are well known present or later-developed, equivalent structures, material or acts for performing the claimed function.

BRIEF DESCRIPTION OF THE DRAWINGS

A more complete understanding of the present invention may be derived by referring to the detailed description when considered in connection with the following illustrative figures. In the figures, like-reference numbers refer to like-elements or acts throughout the figures.

FIG. 1. Block diagram of the IMM for two APAT and SACore models

FIG. 2. Validation region for two tracks and five measurements

FIG. 3. Feasible joint association events for FIG. 2

FIG. 4. APAT Overview

FIG. 5. APAT Model

FIG. 6. APAT Process

FIG. 7. APAT Process

FIG. 8. Radar/Vision—10 Hz/0 Hz CrossRange & Down Range Error

FIG. 9. Radar/Vision—10 Hz/3 Hz CrossRange & Down Range Error

FIG. 10. Radar/Vision—10 Hz/10 Hz CrossRange & Down Range Error

FIG. 11. Radar/Vision—10 Hz/25 Hz CrossRange & Down Range Error

FIG. 12. Performance Improvements Aggregated

FIG. 13. Performance Improvement

FIG. 14. Performance Optimization

FIG. 15.-FIG. 30. Performance Predictions

DETAILED DESCRIPTION

In the following description, and for the purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the various aspects of the invention. It will be understood, however, by those skilled in the relevant arts, that the present invention may be practiced without these specific details. In other instances, known structures and devices are shown or discussed more generally in order to avoid obscuring the invention. In many cases, a description of the operation is sufficient to enable one to implement the various forms of the invention, particularly when the operation is to be implemented in software. It should be noted that there are many different and alternative configurations, devices and technologies to which the disclosed inventions may be applied. The full scope of the inventions is not limited to the examples that are described below.

In one application of the invention, a system and methods are presented to model fully integrated systems in an Algorithm Prototyping and Analysis and Test tool, hereafter referred to as APAT FIGS. 4-7. APAT models include both the deterministic variety, where the output is a point estimate, and the stochastic type; where some measure of uncertainty characterizes the model outputs whose behavior is non-deterministic in that a state does not fully determine its next state.

With respect now to FIG. 7. APAT Process, this figure describes the system level implementation APAT. APAT is divided into four main parts;

• • User Interface (UI) and test setup (10-40; 100, 130, 150, 170);
  • Target Generation (TG);
    • synthetic truth scenarios (60) for the controlled generation of true target trajectories (80); resulting in controlled synthetic sensor reports (110) with no error
    • ground truth scenarios (70) for the controlled generation of ground truth target trajectories (90); resulting in controlled ground truth reports with known error
  • Situation Awareness Core (SACore) (180) data aggregation engine (DA)
  • Performance Analysis (PA) (200) and Monte Carlo Analysis (50)

UI and Test Setup:

Variables that can be set as fusion global parameters as referenced in the text below as <n>:

• • <1-25> Set sensor parameters FIG. 7 (100)
  • <26-40> Set filter noise values for target maneuver FIG. 7 (130)
  • <41-94> Set filter operational parameters FIG. 7 (170)

Begin Variables That Can Be Set... <SetFusionGlobals>
1.  num_sensors = 3;
2.  sensor_zaxis = zeros(1,num_sensors); % Z position of sensor(meters) with respect to center of vehicle
3.  sensor_xaxis = zeros(1,num_sensors);% X position of sensor(meters) with respect to center of vehicle
4.  sensor_xaxis = [0, 0, 0];
5.  sensor_yaxis = zeros(1,num_sensors); % Y position of sensor(meters) with respect to center of vehicle
6.  sensor_yaxis = [0, 0, 0];
7.  sensor_boresight = zeros(1,num_sensors); % Angle of sensor centerline (°) with respect to vehicle X axis
8.  sensor_boresight = [0, 0, 0];
9.  sensor_boresight_corr = zeros(1,num_sensors);% Boresight correction wrt sensor centerline(degrees)
10. range_max(1) = 150;
11. range_max(2) = 150;
12. range_max(3) = 150;
13. angle_max(1) = 20;
14. angle_max(2) = 40;
15. angle_max(3) = 20;
16. total_sensor_tracks(1) = 5;FROM SENSOR
17. total_sensor_tracks(2) = 5;
18. total_sensor_tracks(3) = 5;
19. angle_sigma_sensor_deg(1) = 2.0;COVARIANCE TERMS
20. angle_sigma_sensor_deg(2) = .25;
21. angle_sigma_sensor_deg(3) = 2.0;
22. r_sigma_sensor11 = .005;% range error as a percent of the range for relevant target
23. r_sigma_sensor(1) = .01;% range error as a percent of the range
24. r_sigma_sensor(2) = .10;% range error as a percent of the range
25. r_sigma_sensor(3) = .01;% range error as a percent of the range
26. if filters == 1
27. q(1,1) = .5;
28. q(2,1) = .5;
29. elseif filters == 2
30. q(1,1) = .01;
31. q(1,2) = .5;
32. q(2,1) = .01;
33. q(2,2) = .5;
34. elseif filters == 3
35. q(1,1) = .01;
36. q(1,2) = .5;
37. q(1,3) = .01;
38. q(2,1) = .01;
39. q(2,2) = .01;
40. q(2,3) = .5;
41. deg2rad = pi/180;
42. rad2deg = 180/pi;
43. pass_criteria = 12;
44. filters = 3
45. sensor_measurements = 2;2 dim pos
46. rejectmax(1) = 3;
47. rejectmax(2) = 2;
48. rejectmax(3) = 3;
49. newtrack(1) = 2;
50. newtrack(2) = 2;
51. newtrack(3) = 2;
52. states = 6;2 DIM OUTPUT, POS/VEL/ACCEL
53. sensor_track_length = 14;
54. sensor_raw_number = 15;
55. tent_track_length = 8;
56. firm_tracks = 0;
57. for i = 1:num_sensors
58. firm_tracks = firm_tracks + total_sensor_tracks(i);
59. posx_sigma_init_sensor(1) = 2.0;STARTING UNCERTAINTY INITIALIZING COVARIANCE TERM
60. posy_sigma_init_sensor(1) = 5.0;
61. Rx(1) = posx_sigma_init_sensor(1) * posx_sigma_init_sensor(1);
62. Ry(1) = posy_sigma_init_sensor(1) * posy_sigma_init_sensor(1);
63. velx_sigma_init_sensor(1) = 4.0;
64. vely_sigma_init_sensor(1) = 4.0;
65. Rvx(1) = velx_sigma_init_sensor(1) * velx_sigma_init_sensor(1);
66. Rvy(1) = vely_sigma_init_sensor(1) * vely_sigma_init_sensor(1);
67. accx_sigma_init_sensor(1) = 1.0;
68. accy_sigma_init_sensor(1) = 1.0;
69. Rax(1) = accx_sigma_init_sensor(1) * accx_sigma_init_sensor(1);
70. Ray(1) = accy_sigma_init_sensor(1) * accy_sigma_init_sensor(1);
71. posx_sigma_init_sensor(2) = 5.0;
72. posy_sigma_init_sensor(2) = 1.5;
73. Rx(2) = posx_sigma_init_sensor(2) * posx_sigma_init_sensor(2);
74. Ry(2) = posy_sigma_init_sensor(2) * posy_sigma_init_sensor(2);
75. velx_sigma_init_sensor(2) = 2.0;
76. vely_sigma_init_sensor(2) = 2.0;
77. Rvx(2) = velx_sigma_init_sensor(2) * velx_sigma_init_sensor(2);
78. Rvy(2) = vely_sigma_init_sensor(2) * vely_sigma_init_sensor(2);
79. accx_sigma_init_sensor(2) = 1.0;
80. accy_sigma_init_sensor(2) = 1.0;
81. Rax(2) = accx_sigma_init_sensor(2) * accx_sigma_init_sensor(2);
82. Ray(2) = accy_sigma_init_sensor(2) * accy_sigma_init_sensor(2);
83. posx_sigma_init_sensor(3) = 2.0;
84. posy_sigma_init_sensor(3) = 5.0;
85. Rx(3) = posx_sigma_init_sensor(3) * posx_sigma_init_sensor(3);
86. Ry(3) = posy_sigma_init_sensor(3) * posy_sigma_init_sensor(3);
87. velx_sigma_init_sensor(3) = 4.0;
88. vely_sigma_init_sensor(3) = 4.0;
89. Rvx(3) = velx_sigma_init_sensor(3) * velx_sigma_init_sensor(3);
90. Rvy(3) = vely_sigma_init_sensor(3) * vely_sigma_init_sensor(3);
91. accx_sigma_init_sensor(3) = 1.0;
92. accy_sigma_init_sensor(3) = 1.0;
93. Rax(3) = accx_sigma_init_sensor(3) * accx_sigma_init_sensor(3);
94. Ray(3) = accy_sigma_init_sensor(3) * accy_sigma_init_sensor(3);
.....End of variables that can be set

Now with respect to the variables; sensor parameters above are those values (variables) that represent the desired sensor system to be modeled. They include such parameters as the number of sensors in the system <1>; the relationship of the sensor frames to the vehicle reference plane <2-9> for the subsequent translation and rotation of axes; the expected report limits in terms of down range (meters) and cross range (degrees) <10-15>; report outputs from the sensor interface in terms of the number of tracks expected to be report for each sensor <16-18>; 1-sigma (mean of tested actual population with 1 standard deviation) covariance terms for cross range and down range error and linear error expectations <19-25>. The sensor values set in <1-25> can be adjusted during the modeling process for performance optimization and analysis of the output.

As an example <19> sets the sensor 1 covariance term of angular accuracy as ±2° tested and measure over the range of the sensor, the value can be opened to an expected error of ±3° and run again and compared to the first; the value can be closed to an expected error of ±10 and run again and compared to the other two; or the value can be opened to an expected error of ±3° of sensor 1 and closed on sensor 2 <20> to an expected error of ±0.1° and run again and compared. This ability allows the user to choose and model hardware and sensor pairing scenarios against cost and capability at an integrated systems level, this will be discussed further below in the Performance Analysis section.

To those skilled in the art of targeting it is well known that it is impossible to model the absolute expected error of a maneuvering target, therefore an Interacting Multiple Model is implemented as discussed below. As modeling variables that can be set as inputs to the model run are the target maneuver accelerations willing to be accepted in measurement reports. A two dimensional environment is modeled, therefore anticipate three models, near constant velocity; cross range acceleration; and down range acceleration <26-40>. These q values are in terms of g's maneuver and are introduced into the IMM filter structure as noise values. Further, it is desired to not always model the most complex models, therefore the number of filters wished to be implemented can be set on any pairings of sensors, as an example <26, 29, 34> identify the three possibilities, with subsequent values set as noise for each option chosen. For the discussion below, the latter is selected to be described in detail as it represents the most complex case.

Now with respect to the SetFusionGlobals above and the SACore filter setup; the filter setup is identified as items <41-94>; some items are housekeeping in nature <41,42>; items <43-58> are values set as a function of the model described below. Items <59-94 set the starting uncertainty covariance terms expected from each sensor with respect the 2 dimensional states identified in <52> 6 which defined down range and cross range values for position (2), acceleration (2) and velocity (2). As example, <59, 63, 67> set the expected error of the values of sensor 1 in terms of down range/cross range (x/y) position, velocity and acceleration measurements, the cases are:

• • <59/60 >; <2 meters/5 meters>
  • <63/64>; <4 meters/sec/4 meters/sec>
  • <67/68>; <1 meter/sec²/1 meter/sec²>Sensor (2) and sensor (3) are represented in <71-94> in the same manner.

Target Generation (TG)/Simulation

The APAT/SACore state estimation system is based on efficient recursive filter techniques which estimate the state of a dynamic system from a series of incomplete and noisy measurements. As time progresses, a measurement is made and reported. This measurement must be within the expectation of the next measurement. Therefore, the Kalman filter is based on linear dynamical systems discretized in the time domain. They are modeled on a Markov chain built on linear operators perturbed by Gaussian noise. The state of the system is represented as a vector of real numbers. At each discrete time increment, a linear operator is applied to the state to generate the new state, with some noise mixed in, and optionally some information from the controls on the system if they are known. Then, another linear operator mixed with more noise generates the visible outputs from the hidden state. The Kalman filter may be regarded as analogous to the hidden Markov model, with the key difference that the hidden state variables are continuous (as opposed to being discrete in the hidden Markov model). Additionally, the hidden Markov model can represent an arbitrary distribution for the next value of the state variables, in contrast to the Gaussian noise model that is used for the Kalman filter. There is a strong duality between the equations of the Kalman Filter and those of the hidden Markov model. A review of this and other models is given in Roweis and Ghahramani (1999).

In order to use the Kalman filter to estimate the internal state of a process given only a sequence of noisy observations, one must model the process in accordance with the framework of the Kalman filter.

The scenario blocks (60, 70) represent the setup of controlled scenarios for either synthetic or controlled target measurements. These scenarios and trajectories represent absolute (80) or controlled (90) “truth” with respect to position, velocity and acceleration of the host vehicle and targets observed. Scenarios can include any kinematic scenario such as over taking targets, meeting targets, observing target behavior to include spawning and merging targets, etc. Each of the scenarios can further include conditions such as degraded environment, degraded sensor operation, clutter, skipped reports, or any host of non-linear or un-expected errors that can and do occur in the environment. The sensor output blocks FIG. 7 (110) or (120) are fed into the DA through block (140). Additional discussion for blocks (160, 180, 190) are described in the section below. A second file representing target truth are sent to (200) for the analysis of performance.

SACore and Tracker FIG. 7. (180)

The instrument measurements are passed to the tracking and data fusion algorithms, which will attempt to produce tracks for each of the targets. Each time a track is altered (update, rate aid, track deletion etc.) it sends a track message to the performance analysis routines, which are used to compile statistics on the performance of the tracking algorithms. This is achieved by comparing the track messages with the target's true position, speed etc.

The SACore is the tracking and fusion filter for integrating multi-sensor data to improve the track file of target vehicles. The problem of accurate reporting of measurements on a target is made more difficult because of the presence of clutter and the fact that the target can maneuver with no prior information. The tracking and sensor data fusion filter consists of four main parts:

Interacting Multiple Model (IMM) Estimation Filter:

To estimate the target state, position, velocity and possibly acceleration, a Kalman Filter can be used. The problem is, to accurately estimate the state, the model in the filter must accurately match the maneuvering of the target. Since the movement of the target is not known ahead of time, an accurate model cannot be designed so errors in the state estimation will occur. Adding process noise to model the target maneuvers or using a maneuver detector to adapt the filter has been used in the art, but detection delays and large estimation errors during maneuvers are still a problem. It is generally accepted that the Interacting Multiple Model (IMM) estimator provides superior tracking performance compared to a single Kalman Filter.

The IMM is based on using several models in parallel to estimate the maneuvering target's states. Each Kalman Filter, or another estimation filter, uses a different model for each maneuver, one models a constant velocity target, another models an acceleration in the longitudinal axis while another models an acceleration in the lateral axis. Switching between these models during each sample period is determined probabilistically. Unlike maneuver detection systems where only one filter model is used at a time, the IMM uses all filters. The overall state estimate output is a weighted combination of the estimates from the individual filters. The weighting is based on the likelihood that a filter model is the correct maneuvering target model.

The IMM estimator is a state estimation algorithm that uses Markovian switching coefficients. A system with these coefficients is described by r models, M¹, M², . . . , M^r, and given probabilities of switching between these models. The event that model j (M^j) is in effect during the sampling period ending at time tk, (t_k−1, t_k] will be denoted by M^j(k). The dynamics and measurement for a linear system are given by

x(k)=Φ^j(k,k−1)x(k−1)+G^j(k,k−1)w_j(k−1),  (1)

and

z(k)=H^j(k)x(k)+ν^j(k),  (2)

where x(k) is the system state at time t_k, z(k) is the measurement vector at time t_k, Φ^j(k,k−1) is the state-transition matrix from time t_k−1 to time t_k for M^j(k), G^j(k,k−1) is the noise input matrix, and H^j(k) is the observation matrix for M^j(k). The process noise vector w^j(k−1) and the measurement noise vector v^j(k) are mutually uncorrelated zero-mean white Gaussian processes with covariance matrices Q^j(k−1) and R^j(k) respectively. The initial conditions for the system state under each model j are Gaussian random variables with mean x^j(0) and covariance P^j(0). These prior statistics are assumed known, as also is μ^j(0)=Pr{M^j(0)}, which is the initial probability of model j at t₀.

The model switching is governed by a finite-state Markov chain according to the probability π_ij=Pr{M^j(k)|Mⁱ(k−1)} of switching from Mⁱ(k−1) to M^j(k). The model switching probabilities, π_ij, are assumed known and an example is

$\begin{array}{cc}{\pi }_{\mathrm{ij}}=\left[\begin{array}{cc}\mathrm{.95}& \mathrm{.05}\\ \mathrm{.05}& \mathrm{.95}\end{array}\right].& \left(3\right)\end{array}$

A block diagram of the IMM estimator with only two models, for simplicity, is shown in FIG. 1.

The inputs to the IMM estimator are {circumflex over (x)}¹(k−1|k−1), {circumflex over (x)}²(k−1|k−1), P¹(k−1|k−1), P²(k−1|k−1), and μ^i|j(k−1|k−1), all from the sampling period ending at t_k−1. Where {circumflex over (x)}¹(k−1|k−1) is the state estimate from filter 1 at time t_k−1 using measurements from time t_k−1 and P¹(k−1|k−1) is the corresponding state covariance matrix. Each of the filters use a different mixture of {circumflex over (x)}¹(k−1|k−1) and {circumflex over (x)}²(k−1|k−1) for their input, For r models, this mixing allows the model-conditioned estimates in the current cycle to be computed using r filters rather than r² filters, which greatly decreases the computational burden. The inputs to the filters, {circumflex over (x)}⁰²(k−1|k−1), {circumflex over (x)}⁰²(k−1|k−1), and the corresponding covariance matrices are computed in the Interaction (Mixing) block.

For the filter matched to M^j(k), the inputs are

$\begin{array}{cc}{\hat{x}}^{0j}\left(k-1|k-1\right)=\sum _{i=1}^r{\mu }^{i|j}\left(k-1|k-1\right){\hat{x}}^i\left(k-1|k-1\right)& \left(4\right)\\ P^{0j}\left(k-1|k-1\right)=\sum _{i=1}^r{\mu }^{i|j}\left(k-1|k-1\right)\left\{P^i\left(k-1|k-1\right)+\left[{\hat{x}}^i\left(k-1|k-1\right)-{\hat{x}}^{0j}\left(k-1|k-1\right)\right]*{\left[{\hat{x}}^i\left(k-1|k-1\right)-{\hat{x}}^{0j}\left(k-1|k-1\right)\right]}^T\right\},& \left(5\right)\end{array}$

where the conditional model probability is

$\begin{array}{cc}\begin{array}{c}{\mu }^{i|j}\left(k-1|k-1\right)=\mathrm{Pr}\left\{M^i\left(k-1\right)|M^j\left(k\right),Z_1^{k-1}\right\}\\ =\frac{1}{{\mu }^j\left(k|k-1\right)}{\pi }_{\mathrm{ij}}{\mu }^i\left(k-1|k-1\right),\end{array}& \left(6\right)\end{array}$

and the predicted model probability is

$\begin{array}{cc}{\mu }^j\left(k|k-1\right)=\mathrm{Pr}\left\{M^j\left(k\right)|Z_1^{k-1}\right\}=\sum _{i=1}^r{\pi }_{\mathrm{ij}}{\mu }^i\left(k-1|k-1\right).& \left(7\right)\end{array}$

Using the measurements, z(k), for the filter matched to M^j(k), the updates are computed using the familiar Kalman Filter equations

{circumflex over (x)} ^j(k|k−1)=Φ^j(k,k−1){circumflex over (x)}^0j(k|k−1),  (8)

P ^j(k|k−1)=Φ^j(k,k−1)P^0j(k|k−1)[Φ^j(k,k−1)]^T+G^j(k,k−1)Q^j(k−1)[G^j(k,k−1)]^T,  (9)

v ^j(k)=z(k)−H^j(k){circumflex over (x)}^j(k|k−1),  (10)

S ^j(k)=H^j(k)P^j(k|k−1)[H^j(k)]^T+R^j(k),  (11)

K ^j(k)=P^j(k|k−1)[H^j(k)]^T[S^j(k)]⁻¹,  (12)

{circumflex over (x)} ^j(k|k)={circumflex over (x)}^j(k|k−1)+K^j(k)v^j(k),  (13)

P ^j(k|k)=[I−K^j(k)H^j(k)]P^j(k|k−1),  (14)

where {circumflex over (x)}^j(k|k−1) is the predicted state estimate under M^j(k), P^j(k|k−1) is the corresponding prediction covariance, v^j(k) is the residual, S^j(k) is the residual covariance matrix, K^j(k) is the Kalman gain matrix, {circumflex over (x)}^j(k|k) is the updated state estimate under M^j(k), and P^j(k|k) is the updated covariance matrix.

The likelihood of the filter matched to M^j(k) is defined by Λ^j(k)=f[z(k)|M^j(k),Z₁^k−1], where f[•|•] denotes a conditional density. Using the assumption of Gaussian statistics, the filter residual and the residual covariance, the likelihood is

$\begin{array}{cc}{\Lambda }^j\left(k\right)=\frac{1}{\sqrt{\det \left[2\pi S^j\left(k\right)\right]}}\exp \left\{-{{\frac{1}{2}\left[v^j\left(k\right)\right]}^T\left[S^j\left(k\right)\right]}^{-1}v^j\left(k\right)\right\}.& \left(15\right)\end{array}$

The probability for M^j(k) is

$\begin{array}{cc}{\mu }^j\left(k|k\right)=\mathrm{Pr}\left\{M^j\left(k\right)|Z_1^k\right\}=\frac{1}{c}{\mu }^j\left(k|k-1\right){\Lambda }^j\left(k\right),& \left(16\right)\end{array}$

where the normalization factor c is

$\begin{array}{cc}c=\sum _{j=1}^r{\mu }^i\left(k|k-1\right){\Lambda }^i\left(k\right).& \left(17\right)\end{array}$

These computations are performed in the Model Probability Update block.

Finally the combined state estimate {circumflex over (x)}(k|k) and the corresponding state error covariance for the IMM are given by

$\begin{array}{cc}\hat{x}\left(k|k\right)=\sum _{j=1}^r{\mu }^j\left(k|k\right){\hat{x}}^j\left(k|k\right),& \left(18\right)\\ P\left(k|k\right)=\sum _{j=1}^r{\mu }^j\left(k|k\right)\left\{P^j\left(k|k\right)+{\left[{\hat{x}}^j\left(k|k\right)-\hat{x}\left(k|k\right)\right]\left[{\hat{x}}^j\left(k|k\right)-\hat{x}\left(k|k\right)\right]}^T\right\}.& \left(19\right)\end{array}$

The final state estimate, i(k|k), is the best estimate of the target state and P(k|k) is the error covariance matrix for this optimal state estimate.

Nearest Neighbor Joint Probabilistic Data Association Filter:

This section deals with the problem of taking a target measurement from a sensor and either associating it to an existing track, possibly rejecting the measurement because it is clutter, or sending it to the track management process to start a new track.

The process is to first define a validation region for each track and to identify all sensor measurements that fall within that region. For the two-dimensional case the validation region for a track is constructed around the predicted measurement for the track. Using the previous track state, the predicted measurement is computed by propagating this state to the time of the next measurement. This predicted measurement is the center of the validation region. If the measurement falls within this validation region, then it is considered as a candidate for association to the track; otherwise, it is rejected.

The predicted measurement {circumflex over (z)}(k|k−1) is

{circumflex over (z)}(k|k−1)=H(k){circumflex over (x)}(k|k−1),  (20)

where

{circumflex over (x)}(k|k−1)=Φ(k,k−1){circumflex over (x)}(k−1|k−1).  (21)

The residual, v(k), is the difference between the actual measurement, z(k), and the predicted measurement and is

v(k)=z(k)−{circumflex over (z)}(k|k−1)=z(k)−H(k){circumflex over (x)}(k|k−1).  (22)

The residual covariance S(k) is

S(k)=H(k)P(k|k−1)H^T(k)+R(k).  (23)

The residual statistics are Gaussian with zero-mean and the M-dimensional residual density f(v(k)) is

$\begin{array}{cc}f\left(v\left(k\right)\right)=\frac{1}{\sqrt{\det \left[2\pi S\left(k\right)\right]}}\exp \left\{-\frac{1}{2}d^2\right\},& \left(24\right)\end{array}$

where the normalized (squared) distance d² is defined as

d ² =v ^T(k)S⁻¹(k)v(k).  (25)

The normalized distance d² is a chi-square variant, χ_M², with M degrees of freedom.

A measurement will be in the validation region if

[ρ−{circumflex over (z)}(k|k−1)]^TS⁻¹(k)[ρ−{circumflex over (z)}(k|k−1)]≦g²,  (26)

where g is the threshold that insures that the measurement falls in the validation region with a probability of P_G. The threshold g is referred to as the number of sigmas or standard deviations for the gate, and g is determined from a χ_M² table.

FIG. 2 shows an example of two elliptical validation regions for a system with two tracks and four measurements. T1 and T2 are the centers of validated regions 1 and 2 respectively and represent the predicted measurements for tracks 1 and 2. Measurements M1 and M2 are valid measurements for track 1 while measurements M2 and M3 are valid measurements for track 2. Measurement M4 is not a valid measurement for either track. In both cases there is more than one measurement valid for each track and a procedure is needed to determine which measurements will be used to update the track state.

To save computation time, a coarse gating procedure of testing individual components of the residual vector may be used using rectangular gates. That is, a measurement z(k) falls inside the coarse gate region provided that

|z_i(k)−{circumflex over (z)}_i(k|k−1)≦g√{square root over (S_ii(k))},  (27)

for each i=1, . . . , M, where z_i(k) and {circumflex over (z)}_i(k|k−1) are the i^th components of z(k) and {circumflex over (z)}(k|k−1), respectively, and S_ii(k) is the i^th variance in the residual covariance matrix. The rectangular gating eliminates unlikely measurements-to-track pairing and the ellipsoidal gating is used for the measurements that fall in the rectangular gates.

The data association filter that is used for V2.1 is the nearest-neighbor joint probabilistic data association filter (NNJPDAF). To help explain what the NNJPDAF is, the probabilistic data association filter (PDAF) for the single target case will be discussed first. Next the more complex joint probabilistic data association filter (JPDAF), which is used for multiple targets, will be discussed, and then the nearest-neighbor JPDA.

The PDAF is a Bayesian approach that computes the probability that each measurement in a track's validation region is the correct measurement and the probability that none of the validated measurements is the correct measurement. The non-parametric version of the filter will be used since it does not require prior knowledge of the special density of the clutter.

For m measurements falling inside the validation region at time t_k, the probability that the j^th validated measurement z_j(k) is target originated, denoted β_j, is

$\begin{array}{cc}{\beta }_j=\frac{e_j}{b+{\sum }_{l=1}^me_l}\left(j=1,\dots ,m\right),& \left(28\right)\end{array}$

while the probability that none of the measurements is target originated, denoted by β₀, is

$\begin{array}{cc}{\beta }_0=\frac{b}{b+\sum _{l=1}^me_l}.& \left(29\right)\end{array}$

The term e_j is given by

e _j=exp{−½v_j^T(k)S⁻¹v_j(k)},  (30)

where v_j(k) is the residual for the j^th validated measurement and S(k) is the residual covariance for the measurements. The term b is given by

$\begin{array}{cc}b=\frac{m}{V}\sqrt{\det \left[2\pi S\left(k\right)\right]}\frac{1-P_DP_G}{P_D},& \left(31\right)\end{array}$

where V is the hypervolume of the validated region at t_k, P_D is the detection probability (assumed known), and P_G is the probability of the target originated measurement falling inside the validated region (usually set to at least 0.95). The term b accounts for the possibility that none of the validated measurements is target originated and that the target-originated measurement was not detected, or fell outside of the validated region.

The state in the PDAF is updated using all of the measurements and the combined residual v(k) as follows:

$\begin{array}{cc}\hat{x}\left(kk\right)=\hat{x}\left(kk-1\right)+K\left(k\right)v\left(k\right),& \left(32\right)\\ v\left(k\right)=\sum _{j=1}^m{\beta }_jv_j\left(k\right),& \left(33\right)\end{array}$

and v_j(k) is the residual for the j^th validated measurement,

v _j(k)=z_j(k)−H(k){circumflex over (x)}(k|k−1).  (34)

The updated covariance is given by

P(k|k)=β₀P(k|k−1)+[1−β₀]P^C(k)+{tilde over (P)}(k),  (35)

where

P ^C(k)=P(k|k−1)−K(k)S(k)K^T(k),  (36)

and

$\begin{array}{cc}\tilde{P}\left(k\right)=K\left(k\right)\left[\sum _{j=1}^m{\beta }_jv_j\left(k\right)v_j^T\left(k\right)-v\left(k\right)v^T\left(k\right)\right]K^T\left(k\right).& \left(37\right)\end{array}$

To extend to the multi-target case, the JPDAF is used. First all feasible joint association events θ in the current scan are determined. A feasible event is a set of non-conflicting validated measurement-to-track pairings in which a measurement can originate from only one source, and at most one measurement can originate from a target. Any number of measurements can originate from clutter. Table 1 lists the eight feasible joint association events θ for FIG. 2. The θ_tj's are the single events making up a joint event θ. θ_tj denotes a single event that measurement j (=1, . . . , m) originated from target t (t=1, . . . , N), where m is the total number of measurements in the current scan, N is the total number of targets, and t=0 indicates that the measurement is a clutter detection. In the above example m=3 and N=2. For example, the joint event θ=7 is made up of θ₁₁ (T1 originates M1), θ₀₂ (M2 originates from clutter), and θ₂₃ (T2 originates from M3). The binary target detection indicator δ_t for target t (t=1, . . . , N) has a value of one if a measurement is assigned to target t in θ, and it is zero otherwise. The binary measurement association indicator τ_j for measurement j has a value of one if measurement j is assigned to a target t (t=1, . . . , N) in θ, and is zero otherwise. The quantity φ is the number of measurements originating from clutter in θ.

The joint association event probabilities are given by

$\begin{array}{cc}\mathrm{Pr}\left\{\theta Z_1^k\right\}=\frac{\gamma \left(\theta \right)}{c},& \left(38\right)\end{array}$

where the normalized constant c is

$\begin{array}{cc}c=\sum _{\theta }\gamma \left(\theta \right),& \left(39\right)\end{array}$

where

$\begin{array}{cc}\gamma \left(\theta \right)=\frac{\phi !}{V^{\phi }}\prod _{j=1}^m{\left({\Lambda }_{t_jj}\right)}^{{\tau }_j}\prod _{t=1}^N\left\{{\left(P_D^t\right)}^{{\delta }_t}{\left(1-P_D^t\right)}^{1-{\delta }_t}\right\},& \left(40\right)\end{array}$

and V is the volume of the surveillance region.

The marginal probability of target t (1=0, 1, . . . , N) originating the measurement j, denoted β_ij, is obtained by summing over all feasible joint events θ in which the single event θ_tj occurs, and is given by

$\begin{array}{cc}{\beta }_{\mathrm{tj}}=\mathrm{Pr}\left\{{\theta }_{\mathrm{tj}}Z_1^k\right\}=\sum _{\theta \text{:}{\theta }_{\mathrm{tj}}\in \theta }\mathrm{Pr}\left\{\theta Z_1^k\right\}.& \left(41\right)\end{array}$

Once the marginal association probabilities β_tj are computed, they are used in equations (32) to (37) to update the state for target 1.

For the above example, the normalized constant c is the sum of the eight γ's in Table 1 and is given by

c=P _D ²(Λ₁₁Λ₂₂+Λ₁₁Λ₂₃+Λ₁₂Λ₂₃)V+2!P_D(1−P_D)(Λ₁₁+Λ₁₂+Λ₂₂+Λ₂₃)/V²+3!(1−P_D)2/V³  (42)

The eight joint association probabilities Pr({θ|Z_i^k} are computed by dividing the γ's in Table 1 by c. The marginal association probabilities, β_tj's are then computed using equation (41). As an example, the marginal association probabilities for T1 are

$\begin{array}{cc}\begin{array}{c}{\beta }_{11}=\mathrm{Pr}\left\{{\theta }_{11}Z_1^k\right\}=\mathrm{Pr}\left\{\theta =2Z_1^k\right\}+\\ \mathrm{Pr}\left\{\theta =6Z_1^k\right\}+\mathrm{Pr}\left\{\theta =7Z_1^k\right\}\\ =1/c\left\{2!P_D\left(1-P_D\right){\Lambda }_{11}\right)/V^2+P_D^{2(}{\Lambda }_{11}{\Lambda }_{22}+{\Lambda }_{11}{\Lambda }_{23})V\},\end{array}& \left(43\right)\\ \begin{array}{c}{\beta }_{12}=\mathrm{Pr}\left\{{\theta }_{12}Z_1^k\right\}=\mathrm{Pr}\left\{\theta =3Z_1^k\right\}+\mathrm{Pr}\left\{\theta =8Z_1^k\right\}\\ =1/c\left(2!P_D\left(1-P_D\right){\Lambda }_{12}\right)/V^2+P_D^2{\Lambda }_{12}{\Lambda }_{23}/V\},\end{array}& \left(44\right)\\ {\beta }_{13}=\mathrm{Pr}\left\{{\theta }_{13}Z_1^k\right\}=0,& \left(45\right)\end{array}$

where β₁₃=0 because θ₁₃ does not occur in any of the eight joint events θ; that is, M3 is not validated by T1 in FIG. 2. To obtain the probability β₁₀ of no measurement originating from T1, all joint events in which no measurements are assigned to T1 (i.e., joint events in which δ₁=0) are obtained from Table 1. These joint events are θ=1, 4, 5, so that

β₁₀=Pr{θ=1|Z₁^k}+Pr={θ=4|Z₁^k}+Pr{θ=5|Z₁^k}.  (46)

Since Σ_jβ_tj=Σ_φPr{θ|Z₁^k}=1 for each target t, β₁₀ can be computed more easily using

$\begin{array}{cc}{\beta }_{t0}=1-\sum _{j=1}^m{\beta }_{\mathrm{tj}}.& \left(47\right)\end{array}$

In particular, B_r0 is given by

β₁₀=1−(β₁₁+β₁₂+β₁₃).  (48)

The marginal association probabilities for target T2 can be found in a similar manner.

The technique described above to compute the marginal association probabilities, β_tj's, can only be used when there are few tracks and/or measurements because of the large amount of computation time required. With large numbers of tracks and measurements, a suboptimal method must be used. Two techniques are described in [1]. They are the Cheap JPDA and the Suboptimal JPDA and both techniques, along with the optimal JPDA, are programmed into this filter version. The choice as to which technique to use can be determined at the beginning of the filter run.

The cheap JPDA was developed by Fitzgerald [4]. The probability β_tj of track t associating with measurement j is approximated by

$\begin{array}{cc}{\beta }_{\mathrm{tj}}=\frac{{\Lambda }_{\mathrm{tj}}}{T_t+M_j-{\Lambda }_{\mathrm{tj}}+B}& \left(49\right)\end{array}$

The quantity T₁ is the sum of all likelihoods for target t, which is given by

$\begin{array}{cc}{T}_t=\sum _{j=1}^m{\Lambda }_{\mathrm{tj}}& \left(50\right)\end{array}$

and M_j, which is the sum of all likelihoods for measurement j, is

$M_j=\sum _{t=1}^N{\Lambda }_{\mathrm{tj}}$

The quantity B is a bias to account for clutter and for non-unity probability of detection. Also

$\begin{array}{cc}{\beta }_{t0}=\frac{B}{T_t+B}& \left(52\right)\end{array}$

For the example in FIG. 2, the quantities M_j (j=1, 2, 3) and T_t (t=1,2) are given by,

M ₁=Λ₁₁ T₁=Λ₁₁+Λ₁₂

M ₂=Λ₁₂+Λ₂₂ T₂=Λ₂₂+Λ₂₃

M ₃=Λ₂₃  (53)

The association probabilities for T1 are

$\begin{array}{cc}{\beta }_{10}=\frac{B}{T_1+B}=\frac{B}{{\Lambda }_{11}+{\Lambda }_{12}+B}& \left(54\right)\\ {\beta }_{11}=\frac{{\Lambda }_{11}}{T_1+M_{11}-{\Lambda }_{11}+B}=\frac{{\Lambda }_{11}}{{\Lambda }_{11}+{\Lambda }_{12}+B}& \left(55\right)\\ {\beta }_{12}=\frac{{\Lambda }_{12}}{T_1+M_2-{\Lambda }_{12}+B}=\frac{{\Lambda }_{12}}{{\Lambda }_{11}+{\Lambda }_{12}+{\Lambda }_{22}+B}& \left(56\right)\\ {\beta }_{13}=\frac{{\Lambda }_{13}}{T_1+M_3-{\Lambda }_{13}+B}=0& \left(57\right)\end{array}$

The association probabilities for T2 can be found in a similar manner. Fitzgerald stated that when clutter was significant, a fixed value of B was adequate. Otherwise he found that a value of B=0 gave satisfactory performance.

Although the cheap JPDA is easy to implement and requires much less computational resources, Roecker and Phillis [2] documented some defects in the cheap JPDA. For example, the association probabilities for a given track may not add up to 1. They showed that this may cause a target track to be drawn off by a nearby target track or a nearby false track. To remove some of these defects, the Suboptimal JPDA was developed by Roecker and Phillis [2].

The concept of partial joint events were introduced by Roecker and Phillis to develop the suboptimal JPDA. A partial joint event considers at most two track-to-measurement pairings. That is, a partial joint event consists of(at most) two single events. Also, it assumes the probability of detection for each track is near unity so that all of the targets are detected and all of the joint events have the same number of measurements assigned to clutter. The computational requirement is more than the cheap JPDA, but less than the optimal JPDA.

The association probabilities in the suboptimal JPDA are computed using the following steps:

• • 1. For each track t, form A_t, which is the list of all indices of the validated measurements for track t.
  • 2. For each measurement j, form C_j, which is the list of indices of the tracks which validate measurement j.
  • 3. For each track t, form the union of all track index lists from all of the measurements that are validated by track t while excluding the index of track t. This list of track indices, denoted by L_t, is given by

$\begin{array}{cc}{L}_t=\bigcup _{j\in A_t}C_j-\left\{t\right\}.& \left(58\right)\end{array}$

• • 4. The probability β_tj of track t associating with measurement j is given by

$\begin{array}{cc}{\beta }_{\mathrm{tj}}=\frac{D_{\mathrm{tj}}}{B+{\sum }_{I\in A_t}D_{\mathrm{tl}}},& \left(59\right)\end{array}$

where B is a bias to account for clutter density and

D _tj=Λ_tj if L_t=φ

D _ij=Λ_ijΣ_SεLtN_sj if L_t≠φ  (60)

where φ is the empty set, and

$\begin{array}{cc}{N}_{\mathrm{sj}}=\underset{k\in A_s,k\ne j}{\max }\left\{{\Lambda }_{\mathrm{sk}}\right\}.& \left(61\right)\end{array}$

The probability of no measurements originating from track t is

$\begin{array}{cc}{\beta }_{t0}=\frac{B}{B+{\sum }_{I\in A_t}D_{\mathrm{tl}}}.& \left(62\right)\end{array}$

For the example in FIG. 2, the lists of validated measurements for each track are A₁={1,2} and A₂={2,3}, and the list of tracks validating each measurement are C₁={1}, C₂₁={1,2} and C₃₁={2}. Here T1, M1, etc, are denoted by 1. The track lists L₁ and L₂ are L₁={2} and L₂={1}. The association probabilities for T1 are computed below. A similar procedure can be used to compute the probabilities for T2. The N_sj's needed in these computations are obtained using (61).

$\begin{array}{cc}{N}_{21}=\underset{k\in A_2,k\ne 1}{\max }\left\{{\Lambda }_{2k}\right\}=\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}& \left(63\right)\\ N_{22}=\underset{k\in A_2,k\ne 2}{\max }\left\{{\Lambda }_{2k}\right\}={\Lambda }_{23}& \left(64\right)\\ N_{23}=\underset{k\in A_2,k\ne 3}{\max }\left\{{\Lambda }_{2k}\right\}={\Lambda }_{22}& \left(65\right)\end{array}$

Equation (60) is then used to compute the D_tj's:

$\begin{array}{cc}{D}_{11}={\Lambda }_{11}\sum _{s\in L_1}^{}N_{s1}={\Lambda }_{11}N_{21}={\Lambda }_{11}\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}& \left(66\right)\\ D_{12}={\Lambda }_{12}\sum _{s\in L_1}^{}N_{s2}={\Lambda }_{12}N_{22}={\Lambda }_{12}{\Lambda }_{23}& \left(67\right)\\ D_{13}={\Lambda }_{13}\sum _{s\in L_1}^{}N_{s3}={\Lambda }_{13}N_{23}={\Lambda }_{13}{\Lambda }_{22}=0& \left(68\right)\end{array}$

Finally, (59) and (62) are used to compute the probabilities for T1:

$\begin{array}{cc}\begin{array}{c}{\beta }_{10}=\frac{B}{B+{\sum }_{l\in A_{1t}}D_{1l}}\\ =\frac{B}{{\Lambda }_{11}\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}+{\Lambda }_{12}{\Lambda }_{23}+B}\end{array}& \left(69\right)\\ \begin{array}{c}{\beta }_{11}=\frac{D_{11}}{B+{\sum }_{l\in A_{1t}}D_{1l}}\\ =\frac{{\Lambda }_{11}\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}}{{\Lambda }_{11}\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}+{\Lambda }_{12}{\Lambda }_{23}+B}\end{array}& \left(70\right)\\ \begin{array}{c}{\beta }_{12}=\frac{D_{12}}{B+{\sum }_{l\in A_{1t}}D_{1l}}\\ =\frac{{\Lambda }_{12}{\Lambda }_{23}}{{\Lambda }_{11}\max \left\{{\Lambda }_{22},{\Lambda }_{23}\right\}+{\Lambda }_{12}{\Lambda }_{23}+B}\end{array}& \left(71\right)\\ \begin{array}{c}{\beta }_{13}=\frac{D_{13}}{B+{\sum }_{l\in A_{1t}}D_{1l}}\\ =0.\end{array}& \left(72\right)\end{array}$

To reduce the number of computations required, the nearest-neighbor JPDAF is used instead of the optimal JPDAF. Instead of using equations (32) to (37) to compute the state estimate and covariance matrix, the NNJPDAF looks for the maximum value of β_tj for each measurement-track pair. The corresponding measurement is then used to update the corresponding track. If the maximum β_tj corresponds to track 0, which is clutter, no track is updated. When a track is updated by a measurement, both the track and measurement are no longer used. The procedure is repeated until all associations have been made.

Out-of-Sequence-Measurements (OOSM):

It is very possible that when the measurements are received by the filter from multiple sensors that the measurements do not arrive in the proper time sequence. One reason that this may occur is that the time delay in getting the data from different sensors is not the same. Another reason, and it is probably the major one, is that if the sensor process time, time from obtaining the raw data to the time the sensor data is transmitted, is subtracted from the sensor time, lower frequency sensors will appear out of sequence. Subtracting the process time may be necessary to get the actual time the measurement is valid.

The above techniques are valid for measurements that are in sequence. Different techniques can be used to handle the OOSM. One is to hold the sensor measurements in a buffer and when the OOSM arrives then process it and the stored measurements in the proper sequence. The problem with this technique is that the current track data can get old and is not properly updated until the OOSM arrives. Another technique is to run the filter as if there are no OOSMs which will keep the track data current and then if an OOSM arrives, go back to the OOSM's time and reprocess the filter tracks from this time forward using the OOSM and any measurements after the OOSM. The problem with this technique is that a lot of data must be stored so that the filter can be restarted from the OOSM time and it takes a lot of process time. The technique used in this filter is to run the filter normally with all current data measurements and store only the track covariance matrices at each measurement time. Then when an OOSM arrives, propagate the track state and covariance back in time to the OOSM time, verify its association to a track and then perform the update with the OOSM. Next propagate this updated track state and covariance matrix to the current time.

The procedure is described in FIGS. [5], [6] and [7]. The algorithm to compute the current state and covariance terms using the OOSM is divided into seven steps

• • 1. State Retrodiction (propagating the state backwards in time to the OOSM time)
  • 2. Measurement Retrodiction
  • 3. Mode Likelihood Functions
  • 4. Data Association
  • 5. State Update within each Filter Model
  • 6. Update of the Current Mode Probabilities
  • 7. Update of the Current Combined Estimate and Covariance

State Retrodiction:

The state in each model i of the IMM, i=1, . . . , r, is retrodicted to the time τ of the OOSM, denoted in discrete time as κ. The time τ is l steps lagged, i.e., t_k−l<τ=t_k<t_k−l+1, where t_k is the current time. The retrodiction of the state xⁱ of model i to κ from k is

{circumflex over (x)} ⁱ(κ|k)=Fⁱ(κ|k){circumflex over (x)}ⁱ(k|k)  (73)

where Fⁱ(κ|k) is the backward transition matrix of mode i to κ from k. The covariances associated with the state retrodiction are calculated as follows. The covariance of the cumulative effect of the process noise to k from κ in mode i is

P _vv ⁱ(k,κ|k)=Qⁱ(k,κ)  (74)

For a piecewise constant Wiener process acceleration model, or Discrete Wiener Process Acceleration (DWPA) model, the covariance of the process noise multiplied by the gain matrix Γ is

$\begin{array}{cc}\begin{array}{c}{Q}^i\left(k,\kappa \right)=\Gamma \tilde{q}{\Gamma }^T\\ =\left[\begin{array}{ccc}\frac{1}{4}T^4& \frac{1}{2}T^3& \frac{1}{2}T^2\\ \frac{1}{2}T^3& T^2& T\\ \frac{1}{2}T^2& T& 1\end{array}\right]{\sigma }_v^2\end{array}& \left(75\right)\end{array}$

with T=t_k−t_κ=t_k−τ r and σ_v² the covariance of the process noise.

The cross-covariance between the state at k and the process noise to k from κ in mode i is

P _xv ⁱ(k,κ|k)=Qⁱ(k,κ)−Pⁱ(k|k−1)S*ⁱ(k)⁻¹Qⁱ(k|κ)  (76)

where

S* ⁱ(k)⁻¹=Pⁱ(k|k−1)⁻¹−Pⁱ(k|k−1)⁻¹Pⁱ(k|k)Pⁱ(k|k−1)⁻¹  (77)

Equation (77) is the key to reducing the i-step lag OOSM problem to the 1-step lag problem while keeping the solution nearly optimal. S*ⁱ(k) is the covariance of the “equivalent measurement” at k that summarizes, for model i, all the measurements from k−l+1 to k. To compute this it is necessary to have stored past state covariances Pⁱ(k−l|k−l) for calculation of Pⁱ(k|k−1) for I up to a maximum delay l_max. To compute Pⁱ(k|k−1):

P ⁱ(k|k−1)=Pⁱ(k|k−l)=Fⁱ(k|k−l)Pⁱ(k−l|k−l)Fⁱ(k|k−l)^T+Q(k|k−l)  (78)

where Fⁱ(k|k−l) is the state transition matrix from k−l to k and Q(k|k−l) is the noise contribution as calculated in (75) with T=t_k−t_k−l. The covariance for the state retrodiction is

P ⁱ(κ|k)=Fⁱ(κ|k)[Pⁱ(k|k)+P_vvⁱ(k,κ|k)−P_xvⁱ(k,κ|k)−P_xvⁱ(k,κ|k)^T]Fⁱ(κ,k)^T  (79)

Measurement Retrodiction:

The retrodicted measurement {circumflex over (z)}ⁱ(κ|k) for filter model i of the IMM estimator is calculated based on (73) for each of the r filter models. The retrodicted OOSM for model i is

{circumflex over (z)} ⁱ(κ|k)=H{circumflex over (x)}ⁱ(κ|k)  (80)

The corresponding covariance of the innovation is

S ⁱ(κ|k)=HPⁱ(κ|k)H^T+R(κ)  (81)

Mode Likelihood Functions:

The likelihood function of each mode at time κ is evaluated based on the OOSM and the corresponding retrodiction from the current time. The likelihood function of mode i based on the OOSM from time τ is

κ_i(κ)=|2πSⁱ(κ|k)|^−1/2exp{−½[z(κ)−{circumflex over (z)}ⁱ(κ)]^TSⁱ(κ|k)⁻¹[z(κ)−{circumflex over (z)}ⁱ(κ|k)]}  (82)

The resulting vector consisting of all the likelihood functions of the system modes i=1, . . . , r in the IMM estimator at time κ as

Λ(x)=[Λ_l(κ), . . . ,Λ_r(κ)]^T  (83)

Data Association:

To perform the data association for an OOSM, it is required to have the IMM estimator's overall likelihood function at time κ. For this the mode-conditioned likelihood functions (82) are used, which rely on the mode-conditioned retrodicted measurement and its covariance, and the retrodicted mode probabilities from the current time k to κ. The mode probabilities at time k (before the OOSM) are given by the vector

μ(k|k)=μ(k|Z^k)=[μ_l(k|k), . . . ,μ_r(k|k)]^T  (84)

where Z^k is the cumulative data at k (not including the OOSM z(κ)) and

μ_i(k|k)=Pr{M(k)=i|Z^k}  (85)

Let Π_ij(k₂, k₁) be the Markov chain transition probability matrix between the modes to time t_k2 from time t_k1. The elements of this transition matrix are

Π_ij(k₂,k₁)=Pr{M(k₂)=j|m(k₁)=i}  (86)

And the mode probability vector evolves according to

μ(k₂|k)=Π(k₂,k₁)^Tμ(k₁|k)  (87)

The transition probability matrix to time t_k2 from time t_k1 is

$\begin{array}{cc}\prod \left(k_2,k_1\right)=\frac{1}{{\lambda }_1+{\lambda }_2}\left[\begin{array}{cc}{\lambda }_2+{\lambda }_1{}^{-\left({\lambda }_1+{\lambda }_2\right)T}& {\lambda }_1-{\lambda }_1{}^{-\left({\lambda }_1+{\lambda }_2\right)T}\\ {\lambda }_2-{\lambda }_2{}^{-\left({\lambda }_1+{\lambda }_2\right)T}& {\lambda }_1+{\lambda }_2{}^{-\left({\lambda }_1+{\lambda }_2\right)T}\end{array}\right]& \left(88\right)\end{array}$

where T=|t_k2−t_k1| and

$\frac{1}{{\lambda }_1}\mathrm{and}\frac{1}{{\lambda }_2}$

are the sojourn times in modes 1 and 2. Therefore, the retrodicted mode probability vector to κ from k is

μ(κ|k)=Π(κ|k)^Tμ(k|k)  (89)

The IMM estimator's overall likelihood function for the OOSM is

$\begin{array}{cc}{\Lambda }_{\mathrm{IMM}}\left(\kappa \right)=\sum _{i=1}^r{\Lambda }_i\left(\kappa \right){\mu }_i\left(\kappa k\right)={\Lambda \left(\kappa \right)}^T\mu \left(\kappa k\right)& \left(90\right)\end{array}$

where Λ_i(x) are given in (82) and μ_i(κ|k) follows from (89). This is to be used in the data association to select the delayed measurements to be used in the track update.

State Update within Each Filter Model:

This shows the update of the state at the current time t_k in each model of the IMM using the selected OOSM. The covariance between the state at k and this measurement is calculated as

P _xz ⁱ(k,κ|k)=[Pⁱ(k|k)−P_xvⁱ(k,κ|k)]Fⁱ(κ|k)^TH^T.  (91)

The filter gain used for the update is

W ⁱ(k,κ)=P_xzⁱ(k,κ|k)Sⁱ(κ)⁻¹  (92)

Then update of the most recent state estimate {circumflex over (x)}(k|k) with the OOSM z(κ) is given by

$\begin{array}{cc}\begin{array}{c}{\hat{x}}^i\left(kZ^k,z\left(\kappa \right)\right)={\hat{x}}^i\left(kZ^k\right)\\ ={\hat{x}}^i\left(k\kappa \right)\\ ={\hat{x}}^i\left(k\kappa \right)+W^i\left(k,\kappa \right)\left[z\left(\kappa \right)-{\hat{z}}^i\left(\kappa k\right)\right]\end{array}& \left(93\right)\end{array}$

where the predicted OOSM is given in (80) and the gain is given in (92). The filter-calculated covariance for the updated state estimate is

P ⁱ(k|κ)=Pⁱ(k|k)−P_xzⁱ(k,κ|k)Sⁱ(κ)⁻¹P_xzⁱ(k,κ|k)^T  (94)

Update of the Current Mode Probabilities:

This shows how the mode probabilities are updated at the current time. The conditional mode probabilities at time k (before the OOSM), given by the vector

μ(k|k)=[μ_l(k|k), . . . ,μ_r(k|k)]^T  (95)

are to be updated with the OOSM into

μ(k|κ)=[μ_l(k|κ), . . . ,μ_r(k|κ)]^T  (96)

according to

$\begin{array}{cc}{\mu }_i\left(k\kappa \right)=\frac{1}{c}\left[\sum _{j=1}^r{\Lambda }_k\left(\kappa \right){\Pi }_{\mathrm{ij}}\left(\kappa ,k\right)\right]{\mu }_i\left(kk\right)& \left(97\right)\end{array}$

where the normalized constant is

$\begin{array}{cc}c=\sum _{l=1}^r\sum _{j=1}^r{\Lambda }_j\left(\kappa \right){\Pi }_{\mathrm{ij}}\left(\kappa ,k\right){\mu }_l\left(kk\right).& \left(98\right)\end{array}$

Update of the Current Combined Estimate and Covariance:

Finally the new combined estimate, after the OOSM, is

$\begin{array}{cc}\hat{x}\left(k\kappa \right)=\sum _{j=1}^r{\hat{x}}^j\left(k\kappa \right){\mu }_j\left(k\kappa \right)& \left(99\right)\end{array}$

and its covariance is

P(k|k)=Σμ_j(k|κ){P^j(k|κ)+[{circumflex over (x)}^j(k|κ)−(k|κ)]×[{circumflex over (x)}^j(k|κ)−{circumflex over (x)}(k|κ)]^T}  (100)

Track Management:

A procedure is needed to determine when to initiate a new track and when to drop a track that is no longer being detected by the sensors. Under normal operations, measurements from the sensors are assigned to tracks and are used to update the tracks. If a new measurement is reported by a sensor and it is not associated to an existing track, it may be used to start a new tentative track.

To become a tentative track a sensor must report the same measurement for X number of consecutive scans. The value of X will be determined from real time tests, but initially it is set to 2. If during the next sensor scan the sensor again reports the measurement for the tentative track, the tentative track now becomes a firm track and can be updated by any sensor measurement. If during the next sensor scan the sensor does not report a measurement that corresponds to the tentative track, the tentative track is dropped.

A new firm track is given the lowest integer that has not already been assigned to a firm track. For example, if there are three firm tracks present that have been assigned numbers 1, 2, and 3 then the next firm track will be assigned number 4. If there are three tracks present and they are assigned numbers 1, 4, and 7 then the next firm track will be assigned number 2.

Another function of the Track Management process is to determine if an existing track should be dropped. The Track Management process keeps track of which sensor updates the firm tracks. If all sensors do not update a firm track for at least Y scans, then the firm track is dropped. The value of Y will be determined from real time tests, but initially it is set to 2 for optical sensors and 3 for ranging sensors. As long as one sensor is updating a firm track it will not be dropped.

Sensor Data Conversion:

The sensors usually report their target measurements in a polar reference frame, range and azimuth or bearing. The tracking estimator usually computes its Firm Target File in a Cartesian reference frame because it is easier to understand and visualize. Therefore, the sensor data must be converted to the Cartesian reference frame before it is used by the tracking filter.

The standard conversion from the polar to the Cartesian reference frame is:

x _m =r _m cos θ_m and y_m=r_m sin θ_m,  (101)

where r_m and θ_m are the range and bearing, respectively, of the sensor target in the polar reference frame and x_m and y_m are the downrange and cross range coordinates, respectively, in the converted Cartesian reference frame. However, when dealing with the statistics of the measurements, mean and variance, one cannot use the above equations to transform from the polar to the Cartesian frames. The uncertainty in the range and the bearing is not a perfect ellipsoid so there needs to be debiased correction terms subtracted from (101) to get a better value for the downrange and cross range measurements and their variances.

The following equations give the debiased conversion from a polar coordinate frame to a Cartesian reference frame:

x _dc =r _m cos θ_m−E[{tilde over (x)}|r_m,θ_m],  (102)

y _dc =r _m sin θ_m−E[{tilde over (y)}|r_m,θ_m],  (103)

where

x^dc and y^dc are the final downrange and cross range debiased conversion coordinates of the sensor target and

E[{tilde over (x)}|r _m,θ_m]=r_m cos θ_m(e^−σ−e^−σθ2/2),  (104)

E[{tilde over (y)}|r _m,θ_m]=r_m sin θ_m(e^−σ−e^−σθ2/2).  (105)

The covariance matrix, R_a, for the downrange and cross range coordinates are

R _a ¹¹=var({tilde over (x)}|r_m,θ_m)=r_m²e^−2σθ2[cos² θ_m(cos h 2σ_θ²−cos h σ_θ²)+sin² θ_m(sin h 2σ_θ²−sin h σ_θ²)]+σ_r²e^−2σθ2[cos² θ_m(2 cos h 2σ_θ²−cos h σ_θ²)+sin² θ_m(2 sin h 2σ_θ²−sin σ_θ²)],  (106)

R _a ²²=var({tilde over (y)}|r_m,θ_m)=r_m²e^−2σθ2[sin² θ_m(cos h 2σ_θ²−cos h σ_θ²)+sin² θ_m(sin h 2σ_θ²−sin h σ_θ²)]+σ_r²e^−2σθ2[sin² θ_m(2 cos h 2σ_θ²−cos h σ_θ²)+sin² θ_m(2 sin h 2σ_θ²−sin σ_θ²)],  (107)

R _a ¹² =cov({tilde over (x)},{tilde over (y)},|r_m,θ_m)=sin θ_m cos θ_me^−4σθ2[σ_r²+(r_m²+σ_r²)(1−e^σθ2)],  (108)

where

σ_r² and σ_θ² are the variances of the range and bearing, respectively, in the sensor polar reference frame.

Performance Improvement Analysis

FIG. 15. A sensor that provides range, r, and angle, θ, can be converted to the Cartesian coordinates of down range, x, and cross range, y, by using x=r cos(θ) and y=r sin(θ). However, r=rt+re and θ=θt+θe, where rt stands for true range and re stands for the range error. The same applies to the angle. This means that the down range measurement, xm, and the cross range measurement, ym, are as follows:

xm=xt+xe=(rt+re)cos(θt+θe) and ym=yt+ye=(rt+re)sin(θt+θe).

Expanding these out can give you the error means and covariances, however, it requires knowing the true range, rt, and angle, θt, which are not known. Therefore, the expected values and covariances must be computed based on the measure values of range and angle and not the true values. When this is done the expected value, mean, and covariance of the down range error, xe, and cross range error, ye, give a correction terms for the range and angle terms.

FIG. 16. Assume you have two sensors that take different time intervals to process the raw data. In the slide sensor 1 could be a radar sensor that process the raw data faster and, therefore, has a faster output rate. Sensor 2 could be an optical sensor that takes longer to process the raw data.

Both sensors get the raw data at time t1. Sensor 1 takes the time from t1 to t2 to process its data and then sends it to the filter. Since the process time is known, a time tag can then be applied to the data which was received at t2 with the time t1 which is more accurately the time of the data. For the data received from sensor 1 at t3, the time tag of t2 is given. This is a way to compensate for the sensor latency. Now look at sensor 2 which takes a longer time to process its data than sensor 1. When the data packet arrives at time t2 it is actually given the time stamp of t1 since that more accurately reflects the time of the data. You can see that this will cause a sequence problem if the data rates of the sensors are not the same. Sensor 2's data at its time t2 is tagged as t1 which comes earlier in real time than sensor 1's data at t3 that is tagged as t2. This data is now out of sequence because the filter has processes data from sensor 1 for time t2 and now it gets data from sensor 2 with a time tag of t1. This is where the filter's OOSM (out-of-sequence-measurements) algorithms come into play.

FIG. 17 JPDA is a technique to deal with the problem of taking measurement data and either associating it with an existing track, rejecting it because it may be clutter, or possibly starting a new track. A validation region is generated for each track about the predicted measurement point and all measurements that fall in these regions are considered. The probabilities for all possible combinations of valid measurements to tracks are computed and the highest probabilities are used to match up the measurements and tracks. Presently this approach is used, which is the nearest-neighbor JPDA filter instead of the true optimal JPDA filter which requires more computations.

FIG. 18. Instead of one Kalman Filter to estimate the position, velocity and acceleration of the targets, three filters are proposed, a constant velocity filter, down range acceleration filter, and a cross range acceleration filter. All three filters run at the same time and their outputs are weighted based on what the IMM thinks is the best filter for the trajectory. This allows the constant velocity filter to estimate the target trajectory when the target is not accelerating, the down range acceleration filter to estimate the trajectory when the target is accelerating straight ahead, and the cross range filter to estimate the trajectory when the target is in a turn. Since the best filter is operating during its part of the trajectory, the estimate errors are smaller. Also if the target is accelerating straight ahead and turning, then both the down range and cross range filters may be equally heavily weighted.

FIG. 19. The IMM filters starts off with all three filters (constant velocity, down range acceleration, and cross range acceleration) weighted equally. The IMM switching soon determines that the target is moving at a constant velocity and weights this filter much heavier than the two acceleration filters.

FIG. 20. In this trajectory both vehicles perform a 90 degree constant rate turn to the right. However, the target vehicle is about 25 meters ahead of the host vehicle and one lane to the left so it starts the turn earlier. The start on the plots indicates the starting position.

FIG. 21. These show the relative position, velocity and acceleration of the target vehicle to the host vehicle. Note the acceleration curve which shows the lateral acceleration due to the turning of the vehicles. The spikes occur due to one vehicle turning while the other is not.

FIG. 22. This shows that again initially the three filters start off equally, but when one vehicle turns and the other does not, the filter that looks for cross range acceleration filter is weighted much heavier during the acceleration spike, but decreases gradually after the spike.

FIG. 23. Typical radar values expected

FIG. 24. Typical maneuver scenario

FIG. 25. This shows the relative position of the target vehicle with respect to the host vehicle. The relative velocity is −50 m/sec in the X axis and 0 in the Y axis.

FIG. 26. This shows the relative position and velocity charts for the X (down range) and Y (cross range) axes.

FIG. 27. Using the long range and medium range sensor characteristic data shown in slide 12 raw sensor data is generated. The upper left plot shows the error in the generated down range position for the long range radar and the upper right plots shows the error for the medium range radar. The magenta and yellow lines show the 1-sigma standard deviation values for the noisy data. The maximum range and horizontal FOV determines the start and stop times of the plots.

The bottom plots show the improvement in the position data using the sensor fusion filter. The bottom left plot shows the result of the fusion filter using only position updates. The bottom tight plot uses position and velocity, range-rate, updates. The magenta lines are the true errors since the true track is known. The blue and yellow lines give the filter 1-sigma standard deviation values of the filter error. The decrease in these uncertainty curves for the no velocity update plot is a result of the medium range radar coming into play at about 4.6 seconds. The small bump up after this drop is due to the long range radar dropping out due to FOV. If range-rate updates are available, then the medium range radar does not improve the performance much as shown in the bottom right plot.

Properly tuning the filter should improve the results since the true error points outside the 1-sigma curves are much less than 32%.

FIG. 28. These plots show the filter down range velocity and acceleration error plots when not using velocity updates and when using velocity updates. Again it shows that the filter is not optimally tuned and with tuning the performance should improve slightly.

The increased values of the uncertainty curves when velocity data is used are a result of the long range radar dropping out due to the target reaching its FOV limits.

FIG. 29. This shows the same data as slide 17 except this is for the cross range position instead of the down range position. The decrease in the uncertainty curves is a result of the fact that the cross range errors is equal to the range times the angle uncertainty. As the target approaches the host, the range decreases, and therefore, the cross range error decreases.

The small bump up in the filter—no velocity updates uncertainty curves is a result of the long rang radar dropping out due to the FOV restriction.

FIG. 30. This shows the same as slide 18 except for the cross range values instead of the down range values. The bump up in the bottom right plot is because the inventors did not have the initial covariance values large enough for the acceleration.

The systems described above can use dedicated processor systems, micro controllers, programmable logic devices, or microprocessors that perform some or all of the operations. Some of the operations described above may be implemented in software and other operations may be implemented in hardware.

For the sake of convenience, the operations are described as various interconnected functional blocks or distinct software modules. This is not necessary, however, and there may be cases where these functional blocks or modules are equivalently aggregated into a single logic device, program or operation with unclear boundaries. In any event, the functional blocks and software modules or described features can be implemented by themselves, or in combination with other operations in either hardware or software.

Having described and illustrated the principles of the invention in a preferred embodiment thereof, it should be apparent that the invention may be modified in arrangement and detail without departing from such principles. Claim is made to all modifications and variation coming within the spirit and scope of the following claims.

Claims

1. A sensor system, comprising: a plurality of object detection sensors;a memory;a controller connected to the object detection sensors, wherein the controller is configured to:receive sensor reports from the object detection sensors, wherein the sensor reports include time of detection,determine different mean values from the sensor reports, wherein each of the different mean values is an average value of a central tendency of a one-sigma probability error distribution of a corresponding sensor report in units of cross range and down range values, wherein the one-sigma probability error distribution of a corresponding object detection sensor represents measurement errors and multiple variance values of the detected object,store the sensor reports and the different mean values in memory as stored values,retrieve the stored values from memory,determine a covariance value of a detected object using the retrieved stored values,use a first measured value of the detected object at a first time of detection, the determined covariance value, and a recursive process to make a first estimate of state for the detected object, wherein the first estimate of state includes at least one of position, velocity, and acceleration,store the first estimate of state and determined covariance value in memory,use the recursive process to receive a second measured value of the detected object at a second time of detection,compare the first estimate of state to the second measured value and store the difference in memory.

2. The system of claim 1, wherein the object detection sensors include emitting and non-emitting non-contact sensors.

3. The system of claim 1, wherein data structure of the sensor reports includes latency of the report.

4. The system of claim 1, wherein the object detection sensor reports comprise one or more object attributes of the detected object.

5. The system of claim 4, wherein the one or more object attributes comprise kinematic state of the detected object.

6. The system of claim 1, wherein the structure of the system is a variable structure using interacting multiple models.

7. The system of claim 1, wherein the system uses one of a track oriented or a measurement oriented data association process.

8. The system of claim 1, wherein the memory comprises any hard disk, Read Only Memory (ROM), Dynamic Random Access (RAM) memory, or any combination of different memory devices.

9. The system of claim 1, wherein the controller comprises at least one of a processor, a micro-controller, and a programmable logic device.

10. The system of claim 1, wherein the object detection sensors are mounted to a platform.

11. A state estimation method, comprising: mounting a plurality of object detection sensors to a platform;configuring a controller connected to the object detection sensors to:receive sensor reports from the object detection sensors, wherein the sensor reports include time of detection,determine different mean values from the sensor reports, wherein each of the different mean values is an average value of a central tendency of a one-sigma probability error distribution of a corresponding sensor report in units of cross range and down range values, wherein the one-sigma probability error distribution of a corresponding object detection sensor represents measurement errors and multiple variance values of the detected object,store the sensor reports and the different mean values in memory as stored values,retrieve the stored values from memory,determine a covariance value of a detected object using the retrieved stored values,use a first measured value of the detected object at a first time of detection, the determined covariance value, and a recursive process to make a first estimate of state for the detected object, wherein the first estimate of state includes at least one of position, velocity, and acceleration,store the first estimate of state and determined covariance value in memory,use the recursive process to receive a second measured value of the detected object at a second time of detection,compare the first estimate of state to the second measured value and store the difference in memory.

12. The method of claim 11, wherein the object detection sensors include emitting and non-emitting non-contact sensors.

13. The method of claim 11, wherein data structure of the sensor reports includes latency of the report.

14. The method of claim 11, wherein the object detection sensor reports comprise one or more object attributes of the detected object.

15. The method of claim 14, wherein the one or more object attributes comprise kinematic state of the detected object.

16. The method of claim 11, wherein the structure of the system is a variable structure using interacting multiple models.

17. The method of claim 11, wherein the system uses one of a track oriented or a measurement oriented data association process.

18. The method of claim 11, wherein the memory comprises any hard disk, Read Only Memory (ROM), Dynamic Random Access (RAM) memory, or any combination of different memory devices.

19. The method of claim 11, wherein the controller comprises at least one of a processor, a micro-controller, and a programmable logic device.

20. The method of claim 11, wherein the platform is mounted to a vehicle.