Use Truth Analysis to Track Maneuvering Targets

Since R2025a

In this example, you will learn how to specify a tracker target specification for a maneuvering aircraft based on knowledge gained from analyzing recorded truth. You will then use this tracker target specification to define a multi-object tracker.

The following picture describes the process. The top chevron describes the general process of specifying a task-oriented tracker. The first part of it is specifying what you want to track, which in this case is a maneuvering aircraft. The properties of the maneuvering target specification are listed below. The process on the left follows the Analyze Truth Data and Define Truth Model Example, where beginning with a scenario or recorded truth, you analyze the truth data to extract information regarding the speeds and accelerations of each target. Use these values to specify the maneuvering aircraft spec.

Use Truth Analysis to Specify Tracker Target

The Analyze Truth Data and Define Truth Model Example shows how to analyze truth data, from simulation or recording, in order to determine the target's speed and acceleration limits. In this example, you will use this knowledge to specify a multi-object tracker to track the targets.

The analysis we did let us to arrive at the following conclusions:

Maximum target maneuver is about 70 meters per second squared.
Maximum target speed is about 450 meters per second.
Maximum climb rate is 280 meters per second.

Because the targets perform very sharp maneuvers, you specify the target as a maneuvering aircraft using the trackerTargetSpec function.

tgtspec = trackerTargetSpec("aerospace","aircraft","maneuvering");
disp(tgtspec)

  ManeuveringAircraft with properties:

                 NumManeuvers: 2                
           MaxHorizontalSpeed: [600 600]    m/s 
             MaxVerticalSpeed: [300 300]    m/s 
    MaxHorizontalAcceleration: [10 90]      m/s²
      MaxVerticalAcceleration: [1 90]       m/s²

Furthermore, you specify the speed and acceleration limits on the maneuvering aircraft spec. The specification allows you to define any number of models. In this case, you use two models: The first represents the non-maneuvering high-speed flight regime and the second represents the maneuvering portion of the flight with high accelerations. The speeds for both flight regimes are the same but the accelerations are different.

tgtspec.MaxHorizontalSpeed = [450 450];
tgtspec.MaxVerticalSpeed = [280 280];
tgtspec.MaxHorizontalAcceleration = [10 70];
tgtspec.MaxVerticalAcceleration = [1 70];
disp(tgtspec)

  ManeuveringAircraft with properties:

                 NumManeuvers: 2                
           MaxHorizontalSpeed: [450 450]    m/s 
             MaxVerticalSpeed: [280 280]    m/s 
    MaxHorizontalAcceleration: [10 70]      m/s²
      MaxVerticalAcceleration: [1 70]       m/s²

Define Scenario and Sensor Data

The Analyze Truth Data and Define Truth Model Example defines the scenario, which for the sake of brevity is loaded from a file here.

load("ManeuveringTargetsScenario.mat","scenario");

To visualize the scenario, use the theaterPlot and its plotters as shown in the following block.

thp = theaterPlot(YLimits=[-100 100]*1e3,XLimits=[0 100]*1e3,ZLimits=[-5000 1000]); 
tgp = platformPlotter(thp,DisplayName="Targets",Marker="."); 
rap = platformPlotter(thp,DisplayName="Radar",Marker="o"); 
dep = detectionPlotter(thp,DisplayName="Detections",Marker="x"); 
cvp = coveragePlotter(thp,DisplayName="Coverage"); 
trp = trackPlotter(thp,DisplayName="Tracks",ConnectHistory="on",ColorizeHistory="on");
view(3)

There are six targets in the scenario and one radar platform. Plot their positions and sensor coverage.

poses = platformPoses(scenario);
tgtPoses = vertcat(poses(1:6).Position);
plotPlatform(tgp,tgtPoses);
plotPlatform(rap,poses(7).Position);
plotCoverage(cvp,coverageConfig(scenario));

Figure contains an axes object. The axes object with xlabel X (m), ylabel Y (m) contains 6 objects of type line, patch. One or more of the lines displays its values using only markers These objects represent Targets, Radar, Detections, Coverage, Tracks, (history).

You can run the scenario to simulate the radar data. To save time, you can load the sensor data from a file, which was created using the helperRecordSensorData function attached at the bottom of this script by uncommenting the line of code below.

% sensorData = helperRecordSensorData(scenario,2024);
load("activeRadarData.mat","sensorData");

Specify Tracker and Tracking Metric

Specify the radar sensor specification as an aerospace monostatic radar using the trackerSensorSpec function. The syncSensor2spec helper function attached to this script helps synchronize properties from the fusionRadarSensor object to the radar specification.

radarspec = trackerSensorSpec("aerospace","radar","monostatic");
radarspec = syncSensor2spec(scenario.Platforms{7}.Sensors{1},radarspec);

Define the tracker using the multiSensorTargetTracker function by specifying the target and sensor specifications to be the maneuvering aircraft and monostatic radar specifications, respectively. Due to the high speeds and maneuvers, increase the MaxMahalanobisDistance property to 15.

tracker = multiSensorTargetTracker(tgtspec,radarspec,"jipda");
tracker.MaxMahalanobisDistance = 15;

Finally, you define a tracking metric based on the Optimal Sub-Pattern Assignment (OSPA) metric. You use the OSPA(2) setting to have a comparison of truth trajectories with estimated trajectories over time steps.

ospa2 = trackOSPAMetric("Metric","OSPA(2)","WindowLength",20,"Distance","posabserr","CutoffDistance",500);

Run the Tracker

The following code block allows you to specify the tracker update interval and break the truth and radar data into these intervals.

trackerUpdateInterval = 1.5; % seconds
activeRadarData = helperSplitSensorData(sensorData,scenario.StopTime,trackerUpdateInterval);
pPoses = helperSavePlatformPoses(scenario,trackerUpdateInterval);

The following code resets the tracker, the metric, and the visualization and allows you to run the tracker at different update intervals.

reset(tracker);
reset(ospa2);
clearPlotterData(thp); 
ospaval = zeros(1,numel(activeRadarData));

The main tracking loop is shown in the following code block.

for ind = 1:numel(activeRadarData)
    % Update tracker
    tracks = tracker(activeRadarData(ind));

    % Visualize data and tracks (helper)
    helperPlotActiveRadarData(thp,radarspec,activeRadarData(ind));
    
    [pos,cov] = getTrackPositions(tracks,"constvel");
    plotTrack(trp,pos,cov,string([tracks.TrackID]));
    
    currentPoses = pPoses(1:6,ind); % Exclude the sensor platform.
    plotPlatform(tgp,vertcat(currentPoses.Position));
    drawnow update

    ospaval(ind) = ospa2(tracks,currentPoses);
end

Figure contains an axes object. The axes object with xlabel X (m), ylabel Y (m) contains 15 objects of type line, patch, text. One or more of the lines displays its values using only markers These objects represent Targets, Radar, Detections, Coverage, Tracks, (history).

Analyze Tracking Metric

Observing the visualization above, the tracking quality appears to be good. There are no observable track breaks, swaps, or other issues. However, using a tracking metric is recommended for a quantitative result and if you would like to compare the result with other tracking algorithms or different tracker settings.

figure();
a=axes; 
plot(a,ospaval,'b');
xlabel('Time step');
ylabel('OSPA(2) value');
title('OSPA(2) Track Metric');

Figure contains an axes object. The axes object with title OSPA(2) Track Metric, xlabel Time step, ylabel OSPA(2) value contains an object of type line.

As a reminder, a lower OSPA value indicates better tracking, because OSPA penalizes missing truths, false tracks, redundant tracks, track swaps, and distance between truth and track. The OSPA graph here shows accurate tracking, which was achieved out of the box by specifying the target as maneuvering and setting the appropriate maneuver values.

Summary

This example described the process of configuring a task-oriented tracker using information about the truth maneuvers. The process includes four steps:

Analyze the truth maneuvers.
Extract knowledge about the target speeds and accelerations.
Use that knowledge to configure a tracker target specification.
Specify the tracker with that target specification.

The source of information about the truth can be a scenario or a recording of flight tests. You can use the Tracker Data Importer to import recorded truth and export it as a tracking scenario recording.

Supporting Functions

syncSensor2spec - A helper function to synchronize a sensor model from the scenario to the spec.

function spec = syncSensor2spec(sensor,spec)
spec.MountingLocation = sensor.MountingLocation;
spec.MountingAngles = sensor.MountingAngles;
spec.HasElevation = sensor.HasElevation;
spec.HasRangeRate = sensor.HasRangeRate;
spec.FieldOfView = sensor.FieldOfView;
spec.FieldOfView(1) = 11;
spec.RangeLimits = sensor.RangeLimits;
spec.RangeRateLimits = sensor.RangeRateLimits;
spec.AzimuthResolution = sensor.AzimuthResolution;
spec.RangeResolution = sensor.RangeResolution;
spec.ElevationResolution = sensor.ElevationResolution;
spec.RangeRateResolution = sensor.RangeRateResolution;
spec.DetectionProbability = sensor.DetectionProbability;
spec.FalseAlarmRate = sensor.FalseAlarmRate;
end

helperRecordSensorData - A helper function to record sensor data.

function sensorData = helperRecordSensorData(scenario,seed) %#ok<DEFNU>
if isa(scenario,"trackingScenarioRecording")
    recording = scenario.RecordedData; % We only need the struct
elseif isa(scenario,"struct")
    recording = scenario;
elseif isa(scenario,"trackingScenario") || isa(scenario,"radarScenario")
    disp("Recording the scenario to get sensor data. This may take several minutes.");
    recording = record(scenario,"Rotmat",IncludeSensors=true,RecordingFormat="struct",InitialSeed=seed);
else
    error("This example utility function can only be used with a trackingScenario a trackingScenarioRecording");
end

if ~isfield(recording,"SensorConfigurations") || ~isfield(recording,"Detections")
    error("Scenario recording must contain sensor configurations and detections to record sensor data");
end
% Assume all steps have the same sensors.
sensors = recording(1).SensorConfigurations;
detBuffer = vertcat(recording.Detections);
for i = 1:numel(recording)
    for j = 1:numel(recording(i).SensorConfigurations)
        recording(i).SensorConfigurations(j).Time = recording(i).SimulationTime;
    end
end
configBuffer = reshape([recording.SensorConfigurations],[],1);
sensorData = cell(1,numel(sensors));
for k = 1:numel(sensors)
    sensor = sensors(k);
    sensorIdx = sensor.SensorIndex;
    sensorMountingAngles =  rotmat(recording(1).CoverageConfig(k).Orientation,"frame");
    
    detSensorIndex = cellfun(@(x)x.SensorIndex,detBuffer);
    dets = detBuffer(detSensorIndex == sensorIdx);

    cfgSensorIndex = arrayfun(@(x)x.SensorIndex, configBuffer);
    cfgs = configBuffer(cfgSensorIndex == sensorIdx);
    sensorData{k} = recordMonostaticData(sensorMountingAngles, dets, cfgs);
end
end

recordMonostaticData - A helper function to record sensor data for a monostatic radar.

function activeRadarData = recordMonostaticData(sensorMountingAngles, detections, cfgs)
% Reported config information
lookMountRot = cell(numel(cfgs),1);
for i = 1:numel(cfgs)
    if cfgs(i).MeasurementParameters(1).IsParentToChild
        lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation;
    else
        lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation';
    end
end
beamRot = cellfun(@(x)x*sensorMountingAngles',lookMountRot,UniformOutput=false);
ypr = cellfun(@(x)fusion.internal.frames.rotmat2ypr(x,"degrees"),beamRot,UniformOutput=false);
lookTime = arrayfun(@(x)x.Time,cfgs);
lookAzimuth = cellfun(@(x)x(1),ypr);
lookElevation = cellfun(@(x)-x(2),ypr);

% Reported measurement information
detectionTime = cellfun(@(x)x.Time,detections);
azimuth = cellfun(@(x)x.Measurement(1),detections);
elevation = cellfun(@(x)x.Measurement(2),detections);
range = cellfun(@(x)x.Measurement(3),detections);
rangeRate = cellfun(@(x)x.Measurement(4),detections);
azimuthAccuracy =  cellfun(@(x)x.MeasurementNoise(1,1),detections);
elevationAccuracy =  cellfun(@(x)x.MeasurementNoise(2,2),detections);
rangeAccuracy =  cellfun(@(x)x.MeasurementNoise(3,3),detections);
rangeRateAccuracy =  cellfun(@(x)x.MeasurementNoise(4,4),detections);
platformID =  cellfun(@(x)x.ObjectAttributes{1}.TargetIndex,detections);

falseAlarms = platformID < 0;
toKeep = ~falseAlarms | (rand(numel(falseAlarms),1)<0.01).*falseAlarms;

activeRadarData = struct(LookTime=lookTime(:)',LookAzimuth=lookAzimuth(:)', ...
    LookElevation=lookElevation(:)', DetectionTime=detectionTime(toKeep)', ...
    Azimuth=azimuth(toKeep)',Elevation=elevation(toKeep)', Range=range(toKeep)', ...
    RangeRate=rangeRate(toKeep)',AzimuthAccuracy=sqrt(azimuthAccuracy(toKeep)'), ...
    ElevationAccuracy=sqrt(elevationAccuracy(toKeep)'),RangeAccuracy=sqrt(rangeAccuracy(toKeep)'), ...
    RangeRateAccuracy=sqrt(rangeRateAccuracy(toKeep)'), ...
    PlatformID=platformID(toKeep)');
end

helperSplitSensorData - A helper function to split sensor data between sensor time steps based on the tracker interval.

function varargout=helperSplitSensorData(sensorData,stopTime,trackerInterval)
varargout = cell(1,nargout);
if iscell(sensorData)
    varargout = cellfun(@(s) splitOneSensorData(s,stopTime,trackerInterval), sensorData, UniformOutput=false);
else
    varargout{1} = splitOneSensorData(sensorData,stopTime,trackerInterval);
end
end

splitOneSensorData - A helper to split one sensor data between tracker updates.

function data=splitOneSensorData(sensorData,stopTime,trackerInterval)
frameTime = 0;
ind = 1;
numFrames = ceil(stopTime/trackerInterval);
data = repmat(sensorData,numFrames,1);
while frameTime < stopTime
    withinInterval = sensorData.LookTime >= frameTime & sensorData.LookTime < frameTime+trackerInterval;
    data(ind).LookTime = sensorData.LookTime(withinInterval);
    data(ind).LookAzimuth = sensorData.LookAzimuth(withinInterval);
    data(ind).LookElevation = sensorData.LookElevation(withinInterval);
    withinInterval = sensorData.DetectionTime >=frameTime & sensorData.DetectionTime < frameTime+trackerInterval;
    data(ind).DetectionTime = sensorData.DetectionTime(withinInterval);
    if isfield(sensorData,"Azimuth")
        data(ind).Azimuth = sensorData.Azimuth(1,withinInterval);
        data(ind).AzimuthAccuracy = sensorData.AzimuthAccuracy(1,withinInterval);
    end
    if isfield(sensorData,"Elevation")
        data(ind).Elevation = sensorData.Elevation(1,withinInterval);
        data(ind).ElevationAccuracy = sensorData.ElevationAccuracy(1,withinInterval);
    end
    if isfield(sensorData,"Range")
        data(ind).Range = sensorData.Range(1,withinInterval);
        data(ind).RangeAccuracy = sensorData.RangeAccuracy(1,withinInterval);
    end
    if isfield(sensorData,"RangeRate")
        data(ind).RangeRate = sensorData.RangeRate(1,withinInterval);
        data(ind).RangeRateAccuracy = sensorData.RangeRateAccuracy(1,withinInterval);
    end
    if isfield(sensorData,"PlatformID")
        data(ind).PlatformID = sensorData.PlatformID(1,withinInterval);
    end
    frameTime = frameTime + trackerInterval;
    ind = ind + 1;
end
end

helperSavePlatformPoses - A helper to save platform poses at tracker update times.

function platformPoses = helperSavePlatformPoses(scenario,trackerInterval)
endind = ceil(scenario.StopTime/trackerInterval);
samplePose = struct(PlatformID=1,ClassID=0,Position=zeros(1,3),Velocity=zeros(1,3),Acceleration=zeros(1,3),Orientation=eye(3),AngularVelocity=zeros(1,3));
platformPoses = repmat(samplePose,numel(scenario.Platforms),endind);
sampleTimes = trackerInterval*(1:endind);

for p = 1:numel(scenario.Platforms)
    thisPlatform = scenario.Platforms{p};
    trajectory = thisPlatform.Trajectory;
    if isa(trajectory,"waypointTrajectory") 
        [position,orientation,velocity,acceleration,angularVelocity] = lookupPose(trajectory,sampleTimes);
        orientation = rotmat(orientation,"frame");
        for i = 1:endind
            platformPoses(p,i) = struct(PlatformID=thisPlatform.PlatformID,ClassID=thisPlatform.ClassID, ...
                Position=position(i,:),Velocity=velocity(i,:),Acceleration=acceleration(i,:), ...
                Orientation=orientation(:,:,i),AngularVelocity=angularVelocity(i,:));
        end
    else
        platformPoses(p,:) = repmat(struct(PlatformID=thisPlatform.PlatformID,ClassID=thisPlatform.ClassID, ...
            Position=thisPlatform.Position,Velocity=zeros(1,3),Acceleration=zeros(1,3),Orientation=thisPlatform.Orientation, ...
            AngularVelocity=zeros(1,3)),1,endind);
    end
end
end

helperPlotActiveRadarData - A helper to plot active radar data using the radar spec.

function helperPlotActiveRadarData(thp, radarSpec, radarData)
arguments
    thp (1,1) theaterPlot
    radarSpec
    radarData struct
end

% Find the relevant plotters in the theater plot.
rap = findPlotter(thp,"DisplayName","Radar");
dep = findPlotter(thp,"DisplayName","Detections");
cvp = findPlotter(thp,"DisplayName","Coverage");

% Plot radar platform
plotPlatform(rap, radarSpec.PlatformPosition);

% Compute and plot the coverage configuration
covActive = helperCoverageConfig(radarSpec, radarData);
plotCoverage(cvp, covActive);

% Compute and plot radar measurements
[posMeas,posCov] = radarDataToMeasurementPositions(radarSpec, radarData);
clearData(dep);
plotDetection(dep, posMeas, posCov);
end

radarDataToMeasurementPositions - A helper to compute radar measurement positions.

function [pos,cov] = radarDataToMeasurementPositions(radarSpec, radarData)
n = size(radarData.Azimuth,2);
mp = struct('Frame','spherical',...
    'OriginPosition',zeros(3,1),...
    'OriginVelocity',zeros(3,1),...
    'Orientation',eye(3),...
    'HasAzimuth',true,...
    'HasElevation',true,...
    'HasRange',true,...
    'HasVelocity',true,...
    'IsParentToChild',true);
mp = repmat(mp,2,1);

detection = objectDetection(0,zeros(4,1),'MeasurementParameters',mp,'SensorIndex',1);
pos = zeros(n,3);
cov = zeros(3,3,n);

for i = 1:n
    detection.Measurement = [radarData.Azimuth(i);radarData.Elevation(i);radarData.Range(i);radarData.RangeRate(i)];
    detection.MeasurementNoise = blkdiag(radarData.AzimuthAccuracy(i)^2,radarData.ElevationAccuracy(i)^2,radarData.RangeAccuracy(i)^2,radarData.RangeRateAccuracy(i)^2);
    detection.MeasurementParameters(1).OriginPosition = radarSpec.MountingLocation(:);
    detection.MeasurementParameters(2).OriginPosition = radarSpec.PlatformPosition(:); 
    
    detection.MeasurementParameters(1).OriginVelocity = zeros(3,1);
    detection.MeasurementParameters(2).OriginVelocity = zeros(3,1);

    [~,scanIdx] = min(abs(radarData.LookTime - radarData.DetectionTime(i)));
    lookAng = [radarData.LookAzimuth(scanIdx);radarData.LookElevation(scanIdx)];
    ypr = zeros(1,3,'like',lookAng);
    ypr(1:2) = lookAng;
    ypr(2) = -ypr(2); % Elevation is opposite the sign of pitch
    lookRot = fusion.internal.frames.ypr2rotmat(ypr,'degrees');

    detection.MeasurementParameters(1).Orientation = lookRot*fusion.internal.frames.ypr2rotmat(radarSpec.MountingAngles,'degrees');
    detection.MeasurementParameters(2).Orientation = radarSpec.PlatformOrientation;

    detection.MeasurementParameters(2).Frame = 'rectangular';

    [pos(i,:), ~, cov(:,:,i)] = matlabshared.tracking.internal.fusion.parseDetectionForInitFcn(detection, '', 'double');
end
end

helperCoverageConfig - A helper to compute radar coverage configuration from radar data.

function cvg = helperCoverageConfig(spec, data)
cvg = struct('Index',1,...
    'LookAngle',zeros(2,1),...
    'FieldOfView',zeros(2,1),...
    'ScanLimits',[0 0],...
    'Range',spec.RangeLimits(2),...
    'Position',[0 0 0],...
    'Orientation',quaternion.ones(1,1));
pos = zeros(3,1);
sensorRot = rotmat(quaternion(spec.MountingAngles,'eulerd','ZYX','frame'),'frame');
sensorPos = spec.MountingLocation;
platformRot = spec.PlatformOrientation;
platformPos = spec.PlatformPosition;
pos = fusion.tracker.internal.utils.transformToParent(pos, sensorPos(:), sensorRot);
pos = fusion.tracker.internal.utils.transformToParent(pos, platformPos(:), platformRot);
rot = quaternion(platformRot*sensorRot,'rotmat','frame');
cvg.Position = pos';
cvg.Orientation = rot;
az = data.LookAzimuth;
fovAz = abs(wrapTo180(az(end) - az(1)));
el = data.LookElevation;
fovEl = abs(wrapTo180(el(end) - el(1)));
fov = [fovAz;fovEl] + spec.FieldOfView(:);
idx = max(1,floor(numel(az)/2));
cvg.LookAngle = [wrapTo180(az(idx));wrapTo180(el(idx))];
cvg.FieldOfView = fov;
end