Create an en-route air traffic scenario

You start by creating an Earth-centered scenario.

% Save current random generator state
s = rng;
% Set random seed for predictable results
rng(2020);
% Create scenario
scene = trackingScenario('IsEarthCentered',true,'UpdateRate',1);

Define the airplane model and trajectory

The matfile flightwaypoints attached in this example contains synthesized coordinates and time information of a flight trajectory from Wichita to Chicago. You use a geoTrajectory object to create the flight trajectory.

load('flightwaypoints.mat')
flightRoute = geoTrajectory(lla,timeofarrival);
airplane = platform(scene,'Trajectory',flightRoute);

Nowadays, commercial planes are all equipped with GPS receivers. As the backbone of ADS-B, the accuracy of onboard GPS can be set to conform with ADS-B requirements. The Navigation Accuracy Category used in ADS-B, for position and velocity are referred to as NACp and NACv, respectively. Per FAA regulation[1], the NACp must be less than 0.05 nautical miles and the NACv must be less than 10 meters per second. In this example, you use a gpsSensor model with a position accuracy of 50 m and a velocity accuracy of 10 m/s to configure the adsbTransponder model. You also use a more realistic RCS signature for the airplane, inspired by that of a Boeing 737.

posAccuracy = 50; % meters
velAccuracy = 10; % m/s

gps = gpsSensor('PositionInputFormat','Geodetic','HorizontalPositionAccuracy',posAccuracy,...
    'VerticalPositionAccuracy', posAccuracy,'VelocityAccuracy',10);

adsbTx = adsbTransponder('MW2020','UpdateRate', 1,'GPS', gps,'Category',adsbCategory.Large);

load('737rcs.mat');
airplane.Signatures{1} = boeing737rcs;

Add surveillance stations along the route

There are several models of long-range surveillance radars used by the FAA. The Air Route Surveillance Radar 4 (ARSR-4) is a radar introduced in the 1990s which can provide 3D returns of any 1 square-meter object at a long range of 250 nautical miles (463 kilometers). Most of ARSR-4 radars are located along the borders of the continental United-States, while slightly shorter range radars are mostly located at FAA radar sites on the continent. In this example, a single radar type is modeled following common specifications of an ARSR-4 as listed below:

A platform is added to the scenario for each radar site. The RCS signature of those platforms is set to be $-$50 decibel to avoid creating unwanted radar returns.

By default, the radar detections are reported in the radar mounting platform body frame, which in this case is the local North East Down frame at the position of each radar site. However, in this example you set the DetectionCoordinates property to Scenario in order to output detections in the Earth-Centered Earth-Fixed (ECEF) frame, which allows the tracker to process all the detections from different radar sites in a common frame.

% Model an ARSR-4 radar
updaterate = 1/12;
fov = [360;30.001];
elAccuracy = atan2d(0.9,463); % 900m accuracy @ max range
elBiasFraction = 0.1;

arsr4 = fusionRadarSensor(1,'UpdateRate',updaterate,...
    'FieldOfView',[360 ; 30.001],...
    'HasElevation',true,...
    'ScanMode','Mechanical',...
    'MechanicalAzimuthLimits',[0 360],...
    'MechanicalElevationLimits',[-30 0],...
    'HasINS',true,...
    'HasRangeRate',true,...
    'HasFalseAlarms',false,...
    'ReferenceRange',463000,...
    'ReferenceRCS',0,...
    'RangeLimits', [0 463000],...
    'AzimuthResolution',1.4,...
    'AzimuthBiasFraction',0.176/1.4,...
    'ElevationResolution',elAccuracy/elBiasFraction,...
    'ElevationBiasFraction',elBiasFraction,...
    'RangeResolution', 323,...
    'RangeBiasFraction',116/323,... Accuracy / Resolution
    'RangeRateResolution',100,...
    'DetectionCoordinates','Scenario');

% Add ARSR-4 radars at each ARTCC site
radarsitesLLA = [41.4228  -88.0583  0;...
    40.6989  -89.8258  0;...
    39.2219  -95.2461  0];

for i=1:3
    radar = clone(arsr4);
    radar.SensorIndex = i;
    platform(scene,'Position',radarsitesLLA(i,:),...
        'Signatures',rcsSignature('Pattern',-50),'Sensors',{radar});
end

You use adsbReceiver to model the reception of ADS-B messages. An ADS-B message contains the position measured by the airplane own GPS instrument. The message is usually encoded and broadcasted on 1090 MHz channel for nearby ADS-B receivers to pick up. You define a reception range of 200 km around each surveillance station. In this example you assume that surveillance stations have perfect communication between each other. Therefore, a central receiver picks up broadcasted ADS-B messages within range of at least one station.

% Define central adsbReceiver
adsbRx = adsbReceiver('ReceiverIndex',2);
adsbRange = 2e5;

Visualize the scene

You use the trackingGlobeViewer to display platforms, trajectories, detections, and tracks on Earth.

viewer = trackingGlobeViewer('Basemap','streets-dark',...
    'TrackLabelScale',1.3,'TrackHistoryDepth',4000,...
    'CoverageMode','Coverage');
% Show radar sites
plotPlatform(viewer,scene.Platforms(2:end));
% Show radar coverage
covcon = coverageConfig(scene);
plotCoverage(viewer,covcon,'ECEF','Alpha',0.1);
% Show flight route
plotTrajectory(viewer,flightRoute);

Surveillance radars have an antenna blind cone, sometimes referred to as "cone of silence". It is a volume of space, directly above the radar, that cannot be surveilled due to antenna scanning limitations. Overlapping coverages in networks of radars is a common mitigation strategy for this blind cone region. With an overlapping strategy, however, there can still be areas not fully covered by the network. In this example, the southernmost radar's cone of silence (shown in orange in the picture above) is only partially covered by the adjacent radar in the network. This creates a blind spot where the airplane will not be detected by any radar.

Define a central radar tracker and a track fuser.

Typically, one ARTCC maintains tracks for all objects within its surveillance zone and passes the tracking information to the next ARTCC as the object flies into a new zone. In this example, you define a single centralized tracker for all the radars. You use a trackerGNN object to fuse the radar detections of the plane from multiple radars.

radarGNN = trackerGNN('TrackerIndex',1,...
    'MaxNumTracks',50,...
    'FilterInitializationFcn',@initfilter,...
    'ConfirmationThreshold',[3 5],...
    'DeletionThreshold',[5 5],...
    'AssignmentThreshold',[250 2000]);

You also fuse radar tracks with ADS-B tracks obtained from the ADS-B receiver. To do this, you configure a central trackFuser object. You set the ConfirmationThreshold and DeletionThreshold to account for the update rate difference between the ADS-B receiver rate at 1 Hz, and the radar tracker rate at 1/12 Hz. To allow for at least two radar tracks to be assigned to a central track, the number of hits must then be at least $2\times 12$.

fuser = trackFuser('FuserIndex',3,'MaxNumSources',2,...
    "AssignmentThreshold",[250 500],...
    "StateFusion",'Intersection',...
    'StateFusionParameters','trace',...
    'ConfirmationThreshold',[2 3*12],...
    'DeletionThreshold',[4 4]*12);

Track the flight using radars and ADS-B

In this section, you simulate the scenario and step the radar tracker, track fuser, ADS-B transponder, and receiver.

% Declare some variables
detBuffer = {};
radarTrackLog = {};
fusedTrackLog = {};
adsbTrackLog = {};
images = {};
snapTimes = [84,564,1128, 2083];

% Track labels and colors
adsblabel = "       ADS-B";
radarlabel = "  Radar";
fusedlabel = string(sprintf('%s\n',"","Fused"));
adsbclr = [183 70 255]/255;
radarclr = [255 255 17]/255;
fusedclr = [255 105 41]/255;

while advance(scene)
    time = scene.SimulationTime;

    % Update position of airplane on the globe
    plotPlatform(viewer,airplane,'TrajectoryMode','None');

    % Generate and plot synthetic radar detections
    [dets, configs] = detect(scene);
    dets = postProcessDetection(dets);
    detBuffer = [detBuffer; dets]; %#ok<AGROW>
    plotDetection(viewer,detBuffer,'ECEF');

    % tracks
    adsbTracks = objectTrack.empty;
    radarTracks = objectTrack.empty;
    fusedTracks = objectTrack.empty;

    % Generate ADS-B messages
    truePose = pose(airplane, 'true','CoordinateSystem','Geodetic');
    adsbMessage = adsbTx(truePose.Position, truePose.Velocity);

    % Messages can be received within range of each surveillance station
    if isWithinRange(radarsitesLLA, truePose.Position, adsbRange)
        adsbTracks = adsbRx(adsbMessage, time);
    end

    % Update radar tracker at end of scan
    if all([configs.IsValidTime]) && (~isempty(detBuffer) || isLocked(radarGNN))
        radarTracks = radarGNN(detBuffer,time);
        detBuffer = {};
    end
 
    % Fuse tracks
    if isLocked(fuser) || ~isempty([adsbTracks;radarTracks])
        fusedTracks = fuser([adsbTracks;radarTracks],time);
    end

    % plot tracks only once every radar scan
    if all([configs.IsValidTime])

        labels = [repmat(adsblabel,1,numel(adsbTracks)),...
            repmat(radarlabel,1,numel(radarTracks)),...
            repmat(fusedlabel,1,numel(fusedTracks))];
        
        colors = [repmat(adsbclr, numel(adsbTracks), 1) ;...
            repmat(radarclr, numel(radarTracks), 1);...
            repmat(fusedclr, numel(fusedTracks), 1)];
        
        plotTrack(viewer,[adsbTracks; radarTracks; fusedTracks],'ECEF',...
            'LabelStyle','Custom',"CustomLabel",labels,'Color',colors,...
            'LineWidth',3);
    end

    % Record the estimated airplane data for metrics
    fusedTrackLog = [fusedTrackLog, {fusedTracks}]; %#ok<AGROW>
    radarTrackLog = [radarTrackLog, {radarTracks}]; %#ok<AGROW>
    adsbTrackLog = [adsbTrackLog, {adsbTracks}]; %#ok<AGROW>

    % Move camera and take snapshots
    images = moveCamera(viewer,airplane,time,snapTimes,images);

end

figure
imshow(images{1});

At the beginning of the scenario, the airplane is far from the southernmost surveillance station and no ADS-B message is transmitted. As a result, the airplane is tracked by radar only. Note that ARSR detections are relatively inaccurate in altitude, which is generally acceptable as air traffic controller separate airplanes horizontally rather than vertically. The least accurate detections are generated by radar sites located at longer ranges. These detections are nonetheless associated to the track. In the figure, the white line represents the true trajectory and the yellow line represents the trajectory estimated by the radar tracker. The first leg of the flight is far from any radar and the corresponding detections have high measurement noise. Additionally, the constant velocity motion model does not model the motion well during the initial flight turn after taking off. The radar track is passed to the fuser which outputs a fused track, shown in orange, following closely the radar track. The fused track is different from the radar track because the fuser adds its own process noise to the track state.

figure
imshow(images{2});

In the above snapshot, the airplane is within the ADS-B communication range and a new ADS-B track is established. The fuser processed this new track and improved the accuracy of the fused track.

figure
imshow(images{3});

In the above snapshot, the airplane enters the cone of silence. The radar tracker deletes the track after several updates without any new detections. At this point the fuser relies only on the ADS-B track to estimate the position of the airplane.

figure
imshow(images{4});

As the airplane enters the area covered by the second and third surveillance radar stations, a new radar track is established. Detections from the two radar stations are fused by the radar tracker and the track fuser fuses the new radar track with the ADS-B track.

Analyze results

You compare the recorded track data for radar and ADS-B. The true position and velocity are available from the geoTrajectory. You use the OSPA metric to compare the quality of tracking from radar only, ADS-B only, and from radar fused with ADS-B. You use the default settings of the trackOSPAMetric object to compare NEES (normalized estimation error squared) for track position and detection assignment.

ospa = trackOSPAMetric;
radarospa = zeros(ceil(numel(radarTrackLog)/12),1);

% Compute radar tracker OSPA
for i=1:12:numel(radarTrackLog)
    tracks = radarTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    radarospa(i) = ospa(tracks, truths);
end

% Compute ADS-B OSPA
adsbospa = zeros(numel(adsbTrackLog),1);
reset(ospa);
for i=1:numel(adsbTrackLog)
    tracks = adsbTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    adsbospa(i) = ospa(tracks, truths);
end

% Compute fused OSPA
fusedospa = zeros(numel(fusedTrackLog),1);
reset(ospa);
for i=1:numel(fusedTrackLog)
    tracks = fusedTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    fusedospa(i) = ospa(tracks, truths);
end

% Plot OSPA
figure
hold on
plot((0:12:(numel(radarTrackLog)-1))/60, radarospa(1:12:end), "Color",radarclr);
plot((0:(numel(adsbTrackLog)-1))/60,adsbospa, "Color", adsbclr, 'LineWidth',2);
plot((0:(numel(fusedTrackLog)-1))/60,fusedospa, "Color", fusedclr);
l=legend('Radar','ADS-B','Radar + ADS-B');
l.Color = [0.1 0.1 0.1];
l.TextColor = [ 1 1 1];
xlabel('Time (min)')
ylabel('OSPA')
ax = gca;
grid on;
box on;
ax.Color = [0.1 0.1 0.1];
ax.GridColor = [1 1 1];

The OSPA metric shows improvements obtained from fusing ADS-B tracks with radar tracks. From simulation time 19 min to 25 min, the radar only OSPA is high because the airplane is flying over the blind spot of the radar network. The availability of ADS-B in this area significantly increases the tracking performance, as indicated by the fused OSPA. Additionally, the metric shows two peaks at the beginning, which can be attributed to the poor performance of a constant-velocity filter during the initial turns of the trajectory and the unavailability of ADS-B. Around time 40 min, the ADS-B only OSPA is degraded by the loss of ADS-B availability in the region. In later segments of the simulation, both radar and ADS-B are available. Radar only OSPA is overall worse than ADS-B only. This is because the radars have poor vertical accuracy in comparison to a GPS.

Summary

In this example you have learned how to create an Earth-centered scenario and define trajectories using geodetic coordinates. You also learned how to model Air Route Surveillance Radars and generate synthetic detections. You feed these detections to a multi-object tracker and estimate the position, velocity, and heading of the plane. The tracking performance is improved by adding and fusing of ADS-B information. You modeled ADS-B reports and integrated them to the tracking solution. In this example, only a single flight was modeled. The benefit of ADS-B can be further manifested when modeling multiple flights in scenarios where safe separation distances must be maintained by air traffic controllers.

Supporting functions

initfilter defines the extended Kalman filter used by trackerGNN. Airplane motion is well approximated by a constant velocity motion model. Therefore, use a relatively small process noise to allocate more weight to the dynamics compared to the measurements which are expected to be quite noisy at long ranges.

function filter = initfilter(detection)
filter = initcvekf(detection);
filter.StateCovariance = 100*filter.StateCovariance; % initcvekf uses measurement noise as the default
filter.ProcessNoise = 0.1;
end

postProcessDetection post-processes the detections by two operations:

function detsout = postProcessDetection(detsin)
n = numel(detsin);
keep = zeros(1,n,'logical');
for i=1:n
    meas = detsin{i}.Measurement(1:3)';
    lla = fusion.internal.frames.ecef2lla(meas);
    if lla(3)>0
        keep(i) = true;
    else
        keep(i) = false;
    end
end
detsout = detsin(keep);

velCovTransversal = 100^2;
for i=1:numel(detsout)
    velCov = detsout{i}.MeasurementNoise(4:6,4:6);
    [P,D] = eig(velCov);
    [m,ind] = min(diag(D));
    D = velCovTransversal * eye(3);
    D(ind,ind) = m;
    correctedCov = P*D*P';
    detsout{i}.MeasurementNoise(4:6,4:6) = correctedCov;
end
end

isWithinRange returns true if the airplane position is within ADS-B receiver range of any of the surveillance stations.

function flag = isWithinRange(stationsLLA, pos, range)
cartDistance = fusion.internal.frames.lla2ecef(stationsLLA) - fusion.internal.frames.lla2ecef(pos);
flag = any(vecnorm(cartDistance,2,2) < range);
end

moveCamera specifies new camera positions to follow the airplane and take snapshots.

function images = moveCamera(viewer, airplane, time, snapTimes,images)

if mod(time,12) == 0

    pos = airplane.Position;
    if time == 0
        campos(viewer, pos(1),pos(2), 2e4);
    elseif time == 1100 % Zoom out
        campos(viewer,pos(1),pos(2), 1e5);
    elseif time == 2000 % Zoom in
        campos(viewer,pos(1),pos(2), 2.6e4);
    else % Keep previous height but follow airplane
        campos(viewer,pos(1),pos(2));
    end
end

% Snapshots
if any(time == snapTimes)
    drawnow
    img = snapshot(viewer);
    images = [ images, {img}];
end
end

Reference

[1] Automatic Dependent Surveillance - Broadcast (ADS-B): https://www.faa.gov/about/office_org/headquarters_offices/avs/offices/afx/afs/afs400/afs410/ads-b