multiplatform radar detection fusion -凯发k8网页登录

load a recording of a tracking scenario

the multiplatformradardetectiongeneration mat-file contains a tracking scenario recording previously generated using the following command

recording = record(scene,'includesensors',true,'initialseed',2018,'recordingformat','recording')

where scene is the tracking scenario created in the multiplatform radar detection generation example.

load('multiplatformscenariorecording.mat');

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 1g (g=9.8m/s2) to enable tracking of maneuvering targets in the scenario.

the tracker's assignmentthreshold is set to 50 to enable detections with large range biases (due to atmospheric refraction effects at long detection ranges) to be associated with tracks in the tracker. the deletionthreshold is set to 3 to delete redundant tracks quickly.

enable the hasdetectabletrackidsinput to specify the tracks that are within the field of view of at least one radar since the last update. track logic is only evaluated on 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 a central tracker

the following loop runs the tracking scenario recording until the end of the scenario. for each step forward in the scenario, detections are collected for processing by the central tracker. the tracker is updated with these detections every 2 seconds.

trackupdaterate = 0.5;   % update the tracker every 2 seconds
% create a display to show the true, measured, and tracked positions of the
% detected targets and platforms.
theaterdisplay = helpermultiplatfusiondisplay(recording,'plotassignmentmetrics', true);
% construct an object to analyze assignment and error metrics
tam = trackassignmentmetrics('distancefunctionformat','custom',...
    'assignmentdistancefcn',@truthassignmentdistance,...
    'divergencedistancefcn',@truthassignmentdistance);
% initialize buffers
detbuffer = {};
sensorconfigbuffer = [];
alltracks = [];
detectabletrackids = uint32([]);
assignmenttable = [];
% initialize next tracker update time
nexttrackupdatetime = 2;
while ~isdone(recording)
    % read the next record of the recording.
    [time, trueposes, covcon, dets, senscon, sensplatids] = read(recording);
    
    % buffer all detections and sensor configurations until the next tracker update
    detbuffer = [detbuffer ; dets]; %#ok
    sensorconfigbuffer = [sensorconfigbuffer ; senscon']; %#ok
        
    % follow the trackupdaterate to update the tracker
    if time >= nexttrackupdatetime || isdone(recording)
        
        if isempty(detbuffer)
            lastdetectiontime = time;
        else
            lastdetectiontime = detbuffer{end}.time;
        end
        
        if islocked(tracker)
            % collect list of tracks which fell within at least one radar's field
            % of view since the last tracker update
            predictedtracks = predicttrackstotime(tracker, 'all', lastdetectiontime);
            detectabletrackids = detectabletracks(alltracks, predictedtracks, sensorconfigbuffer);
        end
        
        % update tracker. only run track logic on tracks that fell within at
        % least one radar's field of view since the last tracker update
        [confirmedtracks, ~, alltracks] = tracker(detbuffer, lastdetectiontime, detectabletrackids);
        
        % analyze and retrieve the current track-to-truth assignment metrics
        tam(confirmedtracks, trueposes);
        
        % 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, coverages, and tracks
        theaterdisplay(detbuffer, covcon, confirmedtracks, assignmenttable, trueposes, sensplatids);
        
        % empty buffers
        detbuffer = {};
        sensorconfigbuffer = [];
        
        % update next track update time
        nexttrackupdatetime = nexttrackupdatetime   1/trackupdaterate;
    end
    
end

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 seven rows for seven platforms in the scenario. each track is shown as a horizontal line. 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 is added to the platform's row, when multiple lines are shown at the same time for a single platform, the platform has multiple tracks assigned to it. in these cases, the newer track associated with the platform is considered as redundant.

endtime = assignmenttable.time(end);
assignmenttable(endtime,{'trackid','assignedtruthid','totallength','divergencecount','redundancycount','redundancylength'})

ans=9×6 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      24               7              19                0                  1                  4        
    60 sec      32             nan               7                0                  1                  7        
    60 sec      41               3              10                0                  0                  0        

notice that the platforms which have difficulties in maintaining tracks (platforms 4 and 7) are also the platforms furthest from the radars. this poor tracking performance is attributed to the gaussian distribution assumption for the measurement noise. the assumption 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 covariance ellipses corresponding to actual target distribution and the distribution of the target given by a radar sensor. the sensor is 5 km away from the target with an angular resolution of 5 degrees. the actual measurement uncertainty has a concave shape resulting from the spherical sensor detection coordinate frame in which the radar estimates the target's position.

maxcondnum = 300;
figure;
helperplotlongrangecorrection(maxcondnum)

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. in other words, no eigenvalue in the measurement covariance can be more than 300 times smaller than the covariance's largest eigenvalue. this treatment expands 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.

[confirmedtracks,correctedassignmenttable,ctheaterdisplay] = ...
    runmultiplatfusionsim(recording,tracker,@longrangecorrection);

endtime = correctedassignmenttable.time(end);
correctedassignmenttable(endtime,{'trackid','assignedtruthid','totallength','divergencecount','redundancycount','redundancylength'})

ans=7×6 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      38              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.

alldetections = vertcat(recording.recordeddata.detections);
ctheaterdisplay(alldetections,covcon,confirmedtracks,correctedassignmenttable,trueposes, sensplatids) 
axes(ctheaterdisplay.theaterplot.parent)
legend('off')
xlim([-1000 5000]); ylim([-41000 -36000]); zlim([-5000 0]);
view([-90 90])
axis square
title('jet executing horizontal turn')

zoom in the view in which the jet is 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])

from another view in which the jet is 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([-25000 -9000]); ylim([-31000 -19000]); zlim([-9000 -2000]);
view([-45 10])
title('crossing airliners')

switching the 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, respectively. the use of a very large covariance 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. observing the uncertainty covariance of tracks t07 and t08, you can see that they provide a consistent estimate of platforms p04 and p05, respectively.

xlim([-10000 10000]); ylim([-0 10000]); zlim([-12000 -5000]);
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.

summary

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 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;
end

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 radar sensor 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);
    end
    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)if infov
                    isdetectable(itrk) = infov;
                    k = numsteps; %#ok
                    isensor = numsensors; %#ok
                end
            end
        end
    end
end
trackids = [tracks.trackid]';
trackids = trackids(isdetectable);
end

tracktosensor this function returns the track's position in the sensor's coordinate frame. the track structure is returned by the trackergnn object and the config structure defining the sensor's orientation relative to the track's coordinate frame is returned by the radar 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';
    end
    offset = frames(m).originposition;
    pos = bsxfun(@minus,pos,offset);
    pos = rotmat*pos;
end
end

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;
end
end