Multiplatform Radar Detection Fusion

This example shows how to fuse radar detections from a multiplatform radar network. The network includes two airborne and one ground-based long-range radar platforms. See the Multiplatform Radar Detection Generation example for details. A central tracker processes the detections from all platforms at a fixed update interval. This enables you to evaluate the network's performance against different target types and maneuvers as well as platform configurations and locations.

Load the Simulated Sensor Data and Ground Truth

The MultiplatformRadarDetectionGeneration.mat file contains three variables.

  • truthLog is a struct that contains the ground truth kinematic information of all seven targets and platforms at each step of the simulation (stored in the Truth field). The Configurations field contains orientation, beam position and sensor coordinate frame for all four radar sensors at each step of the simulation.

  • dataLog is a struct that contains the detections reported by each sensor along with the uncertainties and signal to noise ratio (SNR). Detection times correspond to simulation times with at least one detection reported.

  • scene contains the trackingScenario which is used for visualization purposes.

path = fullfile(matlabroot,'examples','fusion','main');

Define Central Tracker

Use the trackerGNN as a central tracker that processes detections received from all radar platforms in the scenario.

The tracker uses the initFilter supporting function to initialize a constant velocity extended Kalman filter for each new track. initFilter modifies the filter returned by initcvekf to match the high target velocities. The filter's process noise is set to 1-G to enable tracking of maneuvering targets in the scenario.

The tracker's AssignmentThreshold is increased from its default value of 30 to 50 to enable detections with large range biases due to atmospheric refraction effects at long detection ranges to be associated by the tracker. The DeletionThreshold is decreased from its default value of 5 to 3 to delete redundant tracks more quickly.

Enable the HasDetectableTrackIDsInput to specify the tracks which fell within at least one radar's field of view since the last update. Track logic is only applied to tracks which had a detection opportunity since the last tracker update.

tracker = trackerGNN('FilterInitializationFcn', @initFilter, ...
        'AssignmentThreshold', 50, 'DeletionThreshold', 3, ...
        'HasDetectableTrackIDsInput', true);

Track Targets by Fusing Detections in Central Tracker

The following loop runs the recorded data until the end of the scenario has been reached. For each step forward in the scenario, detections are collected from the dataLog for processing by the central tracker. The tracker is updated with these detections every 2 seconds.

trackUpdateRate = 1/2;   % Update the tracker every 2 seconds
stop = scene.StopTime;
trackUpdateTimes = (1/trackUpdateRate:1/trackUpdateRate:stop)';

% Values used for tracker update logic.
allTracks = [];

% Create a display to show the true, measured, and tracked positions of the
% detected targets and platforms.
theaterDisplay = helperMultiPlatFusionDisplay(scene, 'XLim', [-45 45], ...
    'YLim', [-45 45], 'PlotAssignmentMetrics', true);

detectableTrackIDs = uint32([]);

% construct an object to analyze assignment and error metrics
tam = trackAssignmentMetrics('DistanceFunctionFormat','custom',...
assignmentTable = [];

for i=1:numel(trackUpdateTimes)
    time = trackUpdateTimes(i);

    % Buffer all detections since last track update
    isDetReady = (dataLog.DetectionTime>time-1/trackUpdateRate) & (dataLog.DetectionTime<=time);
    detBuffer = dataLog.Detections(isDetReady);

    % Get all beam positions since last track update
    configSelect = (truthLog.SimulationTime>time-1/trackUpdateRate) & (truthLog.SimulationTime<=time);
    configsBuffer = truthLog.Configurations(configSelect,:);

    lastDetectionTime = max(dataLog.DetectionTime(isDetReady));
    if isempty(lastDetectionTime)
        lastDetectionTime = time; % There was no detections since last update

    % Collect list of tracks which fell within at least one radar's field
    % of view since last tracker update
    if isLocked(tracker)
        predictedtracks = predictTracksToTime(tracker, 'all', lastDetectionTime);
        detectableTrackIDs = detectableTracks(allTracks,predictedtracks,configsBuffer);

    % Update tracker, only run track logic on tracks that fell within at
    % least one radar's field of view since last tracker update
    [confirmedTracks, ~, allTracks] = tracker(detBuffer,lastDetectionTime,detectableTrackIDs);

    % Analyze and retrieve the current track-to-truth assignment metrics
    truths = truthLog.Truth(find(configSelect,1,'last'),:);

    [~,~] = tam(confirmedTracks, truths);

    % Store assignment metrics in a table
    currentAssignmentTable = trackMetricsTable(tam);
    rowTimes = seconds(time*ones(size(currentAssignmentTable,1),1));
    assignmentTable = cat(1,assignmentTable,table2timetable(currentAssignmentTable,'RowTimes',rowTimes));

    % Update display with detections and tracks


At the end of the scenario, you see that multiple tracks have been dropped and replaced. You can also see the association of tracks to platforms for the duration of the scenario. The plot has 7 rows for each of the 7 platforms in the scenario. Each track is shown as a horizontal line drawn within the row of the platform to which it is assigned. Track numbers are annotated at the beginning of the lines. Whenever a track is deleted, its line stops. Whenever a new track is assigned to a platform, a new line starting at the time when the assigned track is confirmed is added to the platform's row. When multiple lines are shown at the same time for a single platform, this platform has multiple tracks assigned to it. In such cases, the newer track associated with the platform is considered redundant.

endTime = assignmentTable.Time(end);
ans =

  9x6 timetable

     Time     TrackID    AssignedTruthID    TotalLength    DivergenceCount    RedundancyCount    RedundancyLength
    ______    _______    _______________    ___________    _______________    _______________    ________________

    60 sec       1               2              27                0                  0                  0        
    60 sec       7               4              27                0                  0                  0        
    60 sec       8               5              26                0                  0                  0        
    60 sec       9               1              22                0                  0                  0        
    60 sec      11             NaN              10                0                  0                  0        
    60 sec      12               6              20                0                  0                  0        
    60 sec      23               7              20                0                  1                  5        
    60 sec      32             NaN               7                0                  1                  7        
    60 sec      41               3              10                0                  0                  0        

Notice that the platforms which have the most difficulty maintaining tracks (platforms 4 and 7) are also the platforms at the longest ranges from the radars. This poor tracking performance is attributed to the effects of modeling the measurement noise as a Gaussian distribution. This works well for targets at short ranges, but at long ranges, the measurement uncertainty deviates from a Gaussian distribution. The following figure compares the 1-sigma ellipse corresponding to the Gaussian measurement uncertainty of a target at 30 km measured by a radar with a 5 degree angular resolution to the actual measurement uncertainty. The actual measurement uncertainty has a concave shape which results from the spherical coordinate frame in which the radar estimates the target's position.

maxCondNum = 300;

To account for the concave shape of the actual covariance at long ranges, the longRangeCorrection supporting function constrains the condition number of the reported measurement noise. The corrected measurement covariance shown above is constrained to a maximum condition number of 300. No eigenvalue in the measurement covariance can be more than 300 times smaller than the covariance's largest eigenvalue. This has the effect of expanding the measurement noise along the range dimension to better match the concavity of the actual measurement uncertainty.

Simulate with Long-Range Covariance Correction

Rerun the previous simulation using the longRangeCorrection supporting function to correct the reported measurement noise at long ranges.

correctedDdets = longRangeCorrection(dataLog.Detections,maxCondNum);
[confirmedTracks,correctedAssignmentTable,ctheaterDisplay] = ...

endTime = correctedAssignmentTable.Time(end);
ans =

  7x6 timetable

     Time     TrackID    AssignedTruthID    TotalLength    DivergenceCount    RedundancyCount    RedundancyLength
    ______    _______    _______________    ___________    _______________    _______________    ________________

    60 sec       1              2               27                0                  0                  0        
    60 sec       7              4               27                0                  0                  0        
    60 sec       8              5               26                0                  0                  0        
    60 sec       9              1               22                0                  0                  0        
    60 sec      11              7               25                0                  0                  0        
    60 sec      12              6               20                0                  0                  0        
    60 sec      39              3               10                0                  0                  0        

The preceding figure shows that by applying the long-range correction, no track-drops or multiple tracks are generated for the entire scenario. In this case, there is exactly one track for each platform detected by the surveillance network.

xlim([-1 5]); ylim([-41 -36]); zlim([-5 0]);
view([-90 90])
axis square
title('Jet Executing Horizontal Turn')

Zooming in on the jet executing a horizontal turn, the track follows the maneuvering target relatively well, even though the motion model used in this example is constant velocity. Tracking the maneuver could be further improved by using an interacting multiple-model (IMM) filter such as the trackingIMM filter.

view([-60 25])

In another view of the jet executing a horizontal turn, you can see that the altitude is estimated correctly, despite the inaccurate altitude measurements from the sensors. One of the sensors does not report altitude at all, as seen by the large vertical ellipsoids, while the other two sensors underestimate their uncertainty in the altitude.

xlim([-25 -9]); ylim([-31 -19]); zlim([-9 -2]);
view([-45 10])
title('Crossing Airliners')

Switching a point of view to focus on the crossing airliners, the same inaccurate altitude measurements are depicted. Note how the red detections are centered at an altitude of 8 km, while the two airliners fly at altitudes of 3 and 4 km. The use of a very large uncertainty in the altitude allows the tracker to ignore the erroneous altitude reading from the red detections, and keep track of the altitude using the other two radars. Looking at the uncertainty covariance of tracks T07 and T08, you can see that they provide a consistent estimate of platforms P04 and P05, respectively.

xlim([-10 10]); ylim([-0 10]); zlim([-12 -5]);
view([-15 10])
title('Airborne Radar Platforms')

The last plot focuses on the two airborne radar platforms. Each platform is detected by the other platform as well as by the ground radar. The platform trajectories cross each other, separated by 1000 m in altitude, and their tracks are consistent with the ground truth.


This example shows how to process detections collected across multiple airborne and ground-based radar platforms in a central tracker. In this example, you learned how the measurement noise at long ranges is not accurately modeled by a Gaussian distribution. The concavity of the 1-sigma ellipse of the measurement noise at these long ranges results in extremely poor tracking performance with dropped tracks and multiple tracks assigned to a single platform. You also learned how to correct the measurement noise for detections at long ranges to improve the continuity of the reported tracks.

Supporting Functions

initFilter This function modifies the function initcvekf to handle higher velocity targets such as the airliners in the scenario.

function filter = initFilter(detection)
filter = initcvekf(detection);
classToUse = class(filter.StateCovariance);

% Airliners can move at speeds around 900 km/h. The velocity is initialized
% to 0, but will need to be able to quickly adapt to aircraft moving at
% these speeds. Use 900 km/h as 1 standard deviation for the velocity
% noise.
spd = 900*1e3/3600; % m/s
velCov = cast(spd^2,classToUse);
cov = filter.StateCovariance;
cov(2,2) = velCov;
cov(4,4) = velCov;
filter.StateCovariance = cov;

% Set filter's process noise to allow for some horizontal maneuver
scaleAccel = cast(10,classToUse);
Q = blkdiag(scaleAccel^2, scaleAccel^2, 1);
filter.ProcessNoise = Q;

detectableTracks This function returns the IDs for tracks that fell within at least one sensor's field of view. The sensor's field of view and orientation relative to the coordinate frame of the tracks is stored in the array of sensor configuration structs. The configuration structs are returned by the monostaticRadarSensor and can be used to transform track positions and velocities to the sensor's coordinate frame.

function trackIDs = detectableTracks(tracks,predictedtracks,configs)
% Identify which tracks fell within a sensor's field of view

numTrks = size(tracks,1);
[numsteps, numSensors] = size(configs);
allposTrack = zeros(3,numsteps);
isDetectable = false(numTrks,1);
for iTrk = 1:numTrks
    % Interpolate track positions between current position and predicted
    % positions for each simulation step
    posTrack = tracks(iTrk).State(1:2:end);
    posPreditedTrack = predictedtracks(iTrk).State(1:2:end);
    for iPos = 1:3
        allposTrack(iPos,:) = linspace(posTrack(iPos),posPreditedTrack(iPos),numsteps);
    for iSensor = 1:numSensors
        thisConfig = configs(:,iSensor);
        for k = 1:numsteps
            if thisConfig(k).IsValidTime
                pos = trackToSensor(allposTrack(:,k),thisConfig(k));
                % Check if predicted track position is in sensor field of
                % view
                [az,el] = cart2sph(pos(1),pos(2),pos(3));
                az = az*180/pi;
                el = el*180/pi;
                inFov = abs(az)<thisConfig(k).FieldOfView(1)/2 && abs(el) < thisConfig(k).FieldOfView(2)/2;
                if inFov
                    isDetectable(iTrk) = inFov;
                    k = numsteps; %#ok<FXSET>
                    iSensor = numSensors; %#ok<FXSET>


trackIDs = [tracks.TrackID]';
trackIDs = trackIDs(isDetectable);

trackToSensor This function returns the track's position in the sensor's coordinate frame. The track structure is returned by the trackerGNN and the config structure defining the sensor's orientation relative to the track's coordinate frame is returned by the monostaticRadarSensor object.

function pos = trackToSensor(pos,config)

frames = config.MeasurementParameters;
for m = numel(frames):-1:1
    rotmat = frames(m).Orientation;
    if ~frames(m).IsParentToChild
        rotmat = rotmat';
    offset = frames(m).OriginPosition;
    pos = bsxfun(@minus,pos,offset);
    pos = rotmat*pos;

longRangeCorrection This function limits the measurement noise accuracy reported by the radar to not exceed a maximum condition number. The condition number is defined as the ratio of the eigenvalues of the measurement noise to the largest eigenvalue.

When targets are detected at long ranges from a radar, the surface curvature of the uncertainty of the measurement is no longer well approximated by an ellipsoid, but takes on that of a concave ellipsoid. The measurement noise must be increased along the range dimension to account for the concavity, producing a planar ellipse which encompasses the concave ellipsoid. There are several techniques in the literature to address this. Here, the maximum condition number of the measurement noise is limited by increasing the smallest eigenvalues to satisfy the maximum condition number constraint. This has the effect of increasing the uncertainty along the range dimension, producing an ellipse which better encloses the concave uncertainty.

function dets = longRangeCorrection(dets,maxCond)
for m = 1:numel(dets)
    R = dets{m}.MeasurementNoise;
    [Q,D] = eig(R);
    Q = real(Q);
    d = real(diag(D));
    dMax = max(d);
    condNums = dMax./d;
    iFix = condNums>maxCond;
    d(iFix) = dMax/maxCond;
    R = Q*diag(d)*Q';
    dets{m}.MeasurementNoise = R;