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 time, 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 controllers.

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

Generate the true airplane 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 part of their onboard instruments. To model the flight instruments, you set the PoseEstimator property of the airplane platform. 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. 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 an insSensor model with a position accuracy of 100 m and a velocity accuracy of 10 m/s.

onboardINS = insSensor('PositionAccuracy',100,'VelocityAccuracy',10,...
    "RandomStream","mt19937ar with seed");
airplane.PoseEstimator = onboardINS;
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 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 NED 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 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 = monostaticRadarSensor(1,'UpdateRate',updaterate,...
    'FieldOfView',fov,...,
    'HasElevation',true,...
    'ScanMode','Mechanical',...
    'MechanicalScanLimits',[0 360; -30 0],...
    'HasINS',true,...
    'HasRangeRate',true,...
    'HasFalseAlarms',false,...
    'ReferenceRange',463000,...
    'ReferenceRCS',0,...
    '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

Define a central tracker

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 with other sensor information such as ADS-B.

tracker = trackerGNN('FilterInitializationFcn',@initfilter,...
    'ConfirmationThreshold',[3 5],...
    'DeletionThreshold',[5 5],...
    'AssignmentThreshold',[1000 Inf]);

Visualize the scene

You use the helperGlobeViewer object attached in this example to display platforms, trajectories, detections, and tracks on the Earth.

gl = helperGlobeViewer;
setCamera(gl,[28.9176  -95.3388  5.8e5],[0 -30 10]);

% Show radar sites
plotPlatform(gl,scene.Platforms(2:end),'d');
% Show radar coverage
covcon = coverageConfig(scene);
plotCoverage(gl,covcon);
% Show flight route
plotTrajectory(gl,airplane);
% Take a snapshot
snap(gl);

Track the flight using radars only

In this first simulation, you will simulate the scene and track the plane based solely on the long range radars.

useADSB = false;
snapTimes = [300 1600];

% Declare loop variables
detBuffer = {};
tLLA = [];
tHeading = [];
tSpeed = [];
tTime = [];
images = {};

% Set updata rates in seconds
trackupdatetime = 12;
adsbupdatetime  = 1;
arsrupdatetime  = 12;

while advance(scene)
    time = scene.SimulationTime;
    
    % Update position of airplane on the globe
    plotTarget(gl,'airplane',airplane,'^','full');
    
    % Generate radar detections at the defined rate
    if mod(time,arsrupdatetime) == 0 
        % Generate synthetic radar detections
        dets = detect(scene);
        dets = removeBelowGround(dets);
        detBuffer = [detBuffer; dets]; %#ok<AGROW>
    end
    
    % Generate ADS-B detections at the defined rate
    if useADSB && mod(time,adsbupdatetime) == 0
        adsb = adsbDetect(time,airplane);
        detBuffer = [detBuffer; adsb]; %#ok<AGROW>
    end
    
    % Fuse detections in tracker and update tracks on the globe
    if mod(time,trackupdatetime) == 0 && time > 0
        % Update detections on the globe
        plotDetection(gl,detBuffer);
        % Tracker needs detections for first call i.e.
        cond = ~isempty(detBuffer) || ~isLocked(tracker);
        if cond
            tracks = tracker(detBuffer,time);
            if ~isempty(tracks)
                % Record the estimated airplane data
                [tLLA, tSpeed, tHeading, tTime] = ...
                    helperNavigationData(tracks,tLLA, tSpeed, tHeading, tTime);
            end
            plotTrack(gl,tracks);
            detBuffer = {};
        end
    end
    
    % Move camera and take snapshots
    images = moveCamera(gl,time,snapTimes,images);
end

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 representes the true trajectory and the yellow line the trajectory estimated by the tracker.

for i=1:numel(images)
    figure
    imshow(images{i});
end

The first leg of the flight is far from any radar and detections have high measurement noise. Additionally, the constant velocity motion model does not model the motion well during the initial flight turn after take-off. As the airplane enters the coverage area of the second surveillance radar, additional detections correct the track shown on the second snapshot. The track (in yellow) follows the truth (in white) much closely as it enters the coverage area of the second radar.

Track the flight using radars and ADS-B

The utility function adsbDetect uses the estimated position measured by the airplane own navigation instrument modeled by the PoseEstimator property of the platform. The estimated position is usually encoded and broadcasted on 1090 MHz channel for nearby ADS-B receivers to pick up. in the United States, there is almost a full coverage of commercial flights cruising at altitudes of 10,000 meters or lower.

useADSB = true;
snapTimes = [200 1244];
[trackLLA2, trackSpeed2, trackHeading2, trackTime2, images] = ...
    simulateScene(gl,scene,tracker,useADSB,snapTimes);

% Reset random seed
rng(s);

Show snapshots at the same simulation time are shown below as in the previous section.

for i=1:numel(images)
    figure
    imshow(images{i});
end

ADS-B detections are represented in purple color on the globe. The track matches more closely with the ground truth though no significant improvement is noticed while entering the middle radar coverage area.

Results

You compare the recorded track data in both simulations. The true position and velocity are available from the geoTrajectory. In air traffic control, dynamic states of interests are positions, represented by latitude, longitude, and altitude, heading, and groundspeed.

% Query truth at the recorded timestamps for both runs
[trueLLA, ~, trueVel] = lookupPose(flightRoute,tTime);
[trueLLAadsb, ~, trueVeladsb] = lookupPose(flightRoute,trackTime2);

figure
tiledlayout(3,1);
nexttile
hold on
plot(tTime/60,trueLLA(:,3),'LineWidth',2,'LineStyle','--','DisplayName','Truth');
plot(tTime/60,tLLA(:,3),'DisplayName','Surveillance radar only');
plot(trackTime2/60,trackLLA2(:,3),'DisplayName','Surveillance radar + ADS-B');
legend('Location','northeastoutside')
xlabel('Time (minute)');
ylabel('Altitude (meter)')
title('Altitude')

nexttile
hold on
plot(tTime/60,vecnorm(trueVel(:,(1:2)),2,2),'LineWidth',2,'LineStyle','--');
plot(tTime/60,tSpeed);
plot(trackTime2/60,trackSpeed2);
xlabel('Time (minute)');
ylabel('Speed (m/s)')
title('Groundspeed')

nexttile
hold on
plot(tTime/60,atan2d(trueVel(:,2),trueVel(:,1)),'LineWidth',2,'LineStyle','--');
plot(tTime/60,tHeading);
plot(trackTime2/60,trackHeading2);
xlabel('Time (minute)');
ylabel('Heading (degree)')
title('Heading')

The results show that the additional data provided by ADS-B has a significant effect on the altitude estimation as altitude is the least accurate information reported by the radars. Meanwhile, the groundspeed and heading estimate is also improved with the benefit of more frequent updates (every second vs every 12 seconds).

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 used a simple approach to model ADS-B reports and integrate them to tracking solution. In this example, only a single flight was modeled. However, the benefit of ADS-B can be further manifested when modeling multiple flights in scenario where safe separation distances must be maintained by air traffic controllers.

Supporting functions

simulateScene contains the simulation loop which can be run with or without ADS-B. The scenario is updated every second. ADS-B detections are generated every second and ARSR detections are available at the end of every 12 seconds scan. The function outputs the estimated position of the plane as latitude (deg), longitude (deg), and WGS84 altitude (m), the estimated ground speed (m/s), the estimated heading (deg), and the corresponding simulation time (s).

function [tLLA,tSpeed,tHeading,tTime,images] = simulateScene(gl,scene,tracker,useADSB,snapTimes)
% Reset the scenario, globe, and seed
rng(2020);
restart(scene);
airplane = scene.Platforms{1};
release(tracker);
clear(gl);
setCamera(gl,[39.4563 -90.1187 2.356e6]);

% Display initial state of the scene
covcon = coverageConfig(scene);
plotPlatform(gl,scene.Platforms(2:end),'d');
plotCoverage(gl,covcon);

% Declare loop variables
detBuffer = {};
tLLA = [];
tHeading = [];
tSpeed = [];
tTime = [];
images = {};

% Set update rates in seconds
if useADSB
    trackupdatetime = 1;
else
    trackupdatetime = 12;
end
adsbupdatetime  = 1;
arsrupdatetime  = 12;

while advance(scene)
    time = scene.SimulationTime;
    
    % Update position of airplane on the globe
    plotTarget(gl,'airplane',airplane,'^','full');
    
    % Generate radar detections at the defined rate
    if mod(time,arsrupdatetime) == 0 && time > 0
        % Generate synthetic radar detections
        dets = detect(scene);
        dets = removeBelowGround(dets);
        detBuffer = [detBuffer; dets]; %#ok<AGROW>
    end
    
    % Generate ADS-B detections at the defined rate
    if useADSB && mod(time,adsbupdatetime) == 0
        adsb = adsbDetect(time,airplane);
        detBuffer = [detBuffer; adsb]; %#ok<AGROW>
    end
    
    % Update detections on the globe
    plotDetection(gl,detBuffer);
    
    % Fuse detections in tracker and update tracks on the globe
    if mod(time,trackupdatetime) == 0 && time > 0
        % Tracker needs detections for first call
        cond = ~(isempty(detBuffer) && ~isLocked(tracker));
        if cond
            tracks = tracker(detBuffer,time);
            if ~isempty(tracks)
                % Record the estimated airplane data
                [tLLA, tSpeed, tHeading, tTime] = ...
                    helperNavigationData(tracks,tLLA, tSpeed, tHeading, tTime);
            end
            plotTrack(gl,tracks);
            detBuffer = {};
        end
    end
    
    % Move camera and take snapshots
    images = moveCamera(gl,time,snapTimes,images);
    
end
end

adsbDetect generates detection of the airplane based on its on-board sensor suite as modeled by the insSensor.

function detection = adsbDetect(time,airplane,~)
%adsbDetect generate ADS-B detections

estimpose = pose(airplane,'CoordinateSystem','Cartesian');
meas = [estimpose.Position(:) ; estimpose.Velocity(:)];
measnoise = diag([100 100 100 10 10 10].^2);

mparams = struct('Frame','Rectangular',...
    'OriginPositin',zeros(3,1),...
    'OriginVelocity',zeros(3,1),...
    'Orientation',eye(3),... ADS-B does not report orientation
    'IsParentToChild', 1,...
    'HasAzimuth', 1,...
    'HasElevation', 1,...
    'HasRange', 1, ...
    'HasVelocity', 1);

attr ={ struct('TargetIndex',airplane.PlatformID,'SNR',10)};

detection = {objectDetection(time, meas, 'SensorIndex',20,...
    'MeasurementNoise',measnoise,...
    'MeasurementParameters',mparams,...
    'ObjectAttributes',attr)};
end

initfilter defines the extended kalman filter used by trackerGNN. Airplane motion is well approximated by a constant velocity motion model. Therefore a rather small process noise will give 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 = 4*filter.StateCovariance; % initcvekf uses measurement noise as the default
filter.ProcessNoise = 0.02*eye(3);
end

removeBelowGround removes radar detections that lie below the surface of the Earth as these cannot be created by a real radar.

function detsout = removeBelowGround(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);
end

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

function images = moveCamera(gl, time, snapTimes,images)

if time == 120
    setCamera(gl,[37, -97.935, 4.328e4],[0 -23 30]);
end

if time == 1244
    setCamera(gl,[38.8693 -95.6109 1.214e4],[0 -26.8 38.7]);
end

if time == 1445
    setCamera(gl,[37.995, -94.60 6.87e4],[0 -15 2]);
end

if time == 2100
    setCamera(gl, [38.9747 -94.6385 1.3e5], [0 -20 55]);
end

if time == 3000
    setCamera(gl,[41.636 -87.011 6.33e4],[0 -25 270]);
end

% Snaps
if any(time == snapTimes)
    img = snap(gl);
    images = [ images, {img}];
end
end