Simulate and Track En-Route Aircraft in Earth-Centered Scenarios

Open Live Script

This example shows you how to use an Earth-Centered trackingScenario and a geoTrajectory object to model a flight trajectory that spans thousands of kilometers. You use two different models to generate synthetic detections of the airplane: a monostatic radar and ADS-B reports. You use a multi-object tracker to estimate the plane trajectory, compare the tracking performance, and explore the impact that ADS-B provides on the overall tracking quality.

In the United States, the Federal Aviation Administration (FAA) is responsible for the regulation of thousands of flights everyday across the entire national airspace. Commercial flights are usually tracked at all times, from their departure airport to arrival. An air traffic control system is a complex multilevel system. Airport control towers are responsible for monitoring and handling the immediate vicinity of the airport while Air Route Traffic Control Centers (ARTCC) are responsible for long range en-route surveillance across various regions that compose the national airspace.

The capability as well as complexity of air traffic/surveillance radars have increased significantly over the past decades. The addition of transponders on aircraft adds a two-way communication between the radar facility and airplanes which allows for very accurate position estimation and benefits decision making at control centers. In 2020, all airplanes flying above 10,000 feet are required to be equipped with an Automatic Dependent Surveillance Broadcast (ADS-B) transponder to broadcast their on-board estimated position. This message is received and processed by air traffic control centers.

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:

Update rate: 12 sec
Maximum range (1 meter-square target): 463 km
Range Resolution: 323 m
Range Accuracy: 116 m
Azimuth field of view: 360 deg
Azimuth Resolution: 1.4 deg
Azimuth Accuracy: 0.176 deg
Height Accuracy: 900 m

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 segements 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:

It removes radar detections that lie below the surface of the Earth as these cannot be created by a real radar.
It increases the detection noise level for the velocity measurement component normal to the radial component.

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