Air Traffic Control

This example shows how to generate an air traffic control scenario, simulate radar detections from an airport surveillance radar (ASR), and configure a global nearest neighbor (GNN) tracker to track the simulated targets using the radar detections. This enables you to evaluate different target scenarios, radar requirements, and tracker configurations without needing access to costly aircraft or equipment. This example covers the entire synthetic data workflow.

Air Traffic Control Scenario

Simulate an air traffic control (ATC) tower and the motion of the targets in the scenario as platforms. Simulation of the motion of the platforms in the scenario is managed by trackingScenario.

Create a trackingScenario and add the ATC tower to the scenario.

% Create tracking scenario
scenario = trackingScenario;

% Add a stationary platform to model the ATC tower
tower = platform(scenario);

Airport Surveillance Radar

Add an airport surveillance radar (ASR) to the ATC tower. A typical ATC tower has a radar mounted 15 meters above the ground. This radar scans mechanically in azimuth at a fixed rate to provide 360 degree coverage in the vicinity of the ATC tower. Common specifications for an ASR are listed:

Sensitivity: 0 dBsm @ 111 km
Mechanical Scan: Azimuth only
Mechanical Scan Rate: 12.5 RPM
Electronic Scan: None
Field of View: 1.4 deg azimuth, 10 deg elevation
Azimuth Resolution: 1.4 deg
Range Resolution: 135 m

Model the ASR with the above specifications using the monostaticRadarSensor.

rpm = 12.5;
fov = [1.4;10];
scanrate = rpm*360/60;  % deg/s
updaterate = scanrate/fov(1); % Hz

radar = monostaticRadarSensor(1, 'Rotator', ...
    'UpdateRate', updaterate, ...           % Hz
    'FieldOfView', fov, ...                 % [az;el] deg
    'MaxMechanicalScanRate', scanrate, ...  % deg/sec
    'AzimuthResolution', fov(1), ...        % deg
    'ReferenceRange', 111e3, ...            % m
    'ReferenceRCS', 0, ...                  % dBsm
    'RangeResolution', 135, ...             % m
    'HasINS', true, ...
    'DetectionCoordinates', 'Scenario');

% Mount radar at the top of the tower
radar.MountingLocation = [0 0 -15];
tower.Sensors = radar;

Tilt the radar so that it surveys a region beginning at 2 degrees above the horizon. To do this, enable elevation and set the mechanical scan limits to span the radar's elevation field of view beginning at 2 degrees above the horizon. Because trackingScenario uses a North-East-Down (NED) coordinate frame, negative elevations correspond to points above the horizon.

% Enable elevation scanning
radar.HasElevation = true;

% Set mechanical elevation scan to begin at 2 degrees above the horizon
elFov = fov(2);
tilt = 2; % deg
radar.MechanicalScanLimits(2,:) = [-fov(2) 0]-tilt; % deg

Set the elevation field of view to be slightly larger than the elevation spanned by the scan limits. This prevents raster scanning in elevation and tilts the radar to point in the middle of the elevation scan limits.

radar.FieldOfView(2) = elFov+1e-3;

The monostaticRadarSensor models range and elevation bias due to atmospheric refraction. These biases become more pronounced at lower altitudes and for targets at long ranges. Because the index of refraction changes (decreases) with altitude, the radar signals propagate along a curved path. This results in the radar observing targets at altitudes which are higher than their true altitude and at ranges beyond their line-of-sight range.

Add three airliners within the ATC control sector. One airliner approaches the ATC from a long range, another departs, and the third is flying tangential to the tower. Model the motion of these airliners over a 60 second interval.

trackingScenario uses a North-East-Down (NED) coordinate frame. When defining the waypoints for the airliners below, the z-coordinate corresponds to down, so heights above the ground are set to negative values.

% Duration of scenario
sceneDuration = 60; % s

% Inbound airliner
ht = 3e3;
spd = 900*1e3/3600; % m/s
wp = [-5e3 -40e3 -ht;-5e3 -40e3+spd*sceneDuration -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Outbound airliner
ht = 4e3;
spd = 700*1e3/3600; % m/s
wp = [20e3 10e3 -ht;20e3+spd*sceneDuration 10e3 -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Tangential airliner
ht = 4e3;
spd = 300*1e3/3600; % m/s
wp = [-20e3 -spd*sceneDuration/2 -ht;-20e3 spd*sceneDuration/2 -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Create a display to show the true, measured, and tracked positions of the
% airliners.
[theater,fig] = helperATCExample('Create Display',scenario,tower,radar);
helperATCExample('Update Display',theater,scenario,tower,radar);

GNN Tracker

Create a trackerGNN to form tracks from the radar detections generated from the three airliners. Update the tracker with the detections generated after the completion of a full 360 degree scan in azimuth.

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 target velocities and tracker update interval.

tracker = trackerGNN( ...
    'Assignment', 'Auction', ...
    'AssignmentThreshold',50, ...
    'FilterInitializationFcn',@initFilter);

Simulate and Track Airliners

The following loop advances the platform positions until the end of the scenario has been reached. For each step forward in the scenario, the radar generates detections from targets in its field of view. The tracker is updated with these detections after the radar has completed a 360 degree scan in azimuth.

% Set simulation to advance at the update rate of the radar
scenario.UpdateRate = radar.UpdateRate;

% Create a buffer to collect the detections from a full scan of the radar
scanBuffer = {};

% Initialize the track array
tracks = [];

% Set random seed for repeatable results
rng(2018)

while advance(scenario) && ishghandle(fig)

    % Current simulation time
    simTime = scenario.SimulationTime;

    % Target poses in the ATC's coordinate frame
    targets = targetPoses(tower);

    % Use the tower's true position as its INS measurement
    ins = pose(tower, 'true');

    % Generate detections on target's in the radar's current field of view
    [dets,~,config] = radar(targets,ins,simTime);

    scanBuffer = [scanBuffer;dets]; %#ok<AGROW>

    % Update tracks when a 360 degree scan is complete
    if config.IsScanDone
        % Update tracker
        tracks = tracker(scanBuffer,simTime);

        % Clear scan buffer for next scan
        scanBuffer = {};
    elseif isLocked(tracker)
        % Predict tracks to the current simulation time
        tracks = predictTracksToTime(tracker,'confirmed',simTime);
    end

    % Update display with current beam position, buffered detections, and
    % track positions
    helperATCExample('Take Snapshots',fig,scenario,config);
    helperATCExample('Update Display',theater,scenario,tower,radar,scanBuffer,tracks);
    helperATCExample('Take Snapshots',fig,scenario,config);

end
helperATCExample('Take Snapshots',fig,scenario,config);

Use helperATCExample to show the first snapshot taken at the completion of the radar's second scan.

helperATCExample('Show Snapshots',1);

The preceding figure shows the scenario at the end of the radar's second 360 degree scan. Radar detections, shown as blue dots, are present for each of the simulated airliners. At this point, the tracker has already been updated by one complete scan of the radar. Internally, the tracker has initialized tracks for each of the airliners. These tracks will be shown after the update following this scan, when the tracks are promoted to confirmed, meeting the tracker's confirmation requirement of 2 hits out of 3 updates.

The next two snapshots show tracking of the outbound airliner.

helperATCExample('Show Snapshots',[2 3]);

The previous figures show the track picture before (on the left) and immediately after (on the right) the tracker updates after the radar's second scan. The detection in the figure on the left is used to update and confirm the initialized track from the previous scan's detection for this airliner. The figure on the right shows the confirmed track's position and velocity. The estimated track velocity is indicated by the black line extending from the track's location. The uncertainty of the track's position estimate is shown as the gray ellipse. After only two detections, the tracker has established an accurate estimate of the outbound airliner's position and velocity. The airliner's true altitude is 4 km and it is traveling east at 700 km/hr.

helperATCExample('Show Snapshots',[4 5]);

The state of the outbound airliner's track is coasted to the end of the third scan and shown in the figure above on the left along with the most recent detection for the airliner. Notice how the track's uncertainty has grown since it was updated in the previous figure. The track after it has been updated with the detection is shown in the figure on the right. You notice that the uncertainty of the track position is reduced after the update. The track uncertainty grows between updates and is reduced whenever it is updated with a new measurement. You also observe that after the third update, the track lies on top of the airliner's true position.

helperATCExample('Show Snapshots',6);

The final figure shows the state of all three airliners' tracks at the end of the scenario. There is exactly one track for each of the three airliners. The same track numbers are assigned to each of the airliners for the entire duration of the scenario, indicating that none of these tracks were dropped during the scenario. The estimated tracks closely match the true position and velocity of the airliners.

truthTrackTable = tabulateData(scenario, tracks) %#ok<NOPTS>

truthTrackTable =

  3x4 table

    TrackID        Altitude              Heading               Speed      
               True    Estimated    True    Estimated    True    Estimated
    _______    _________________    _________________    _________________

     "T1"      4000      4119        90         90       700        701   
     "T2"      4000      4029         0        359       300        285   
     "T3"      3000      3082         0        359       900        890

Visualize tracks in 3D to get a better sense of the estimated altitudes.

axis square
view(50,10)
xlim([-25 35])
ylim([-40 20])
zlim([-6 0])

Visualize Tracks on a Map

Plot the latitude-longitude of each track position in a geographic axes using geoplot. You can use the ned2geodetic function from Mapping Toolbox™ to transform track positions from NED Cartesian coordinates to geodetic coordinates. The latitude and longitude of all tracks are stored in the ATCExample.mat file.

helperATCExample('geoPlotTracks','ATCExample.mat');

Summary

This example shows how to generate an air traffic control scenario, simulate radar detections from an airport surveillance radar (ASR), and configure a global nearest neighbor (GNN) tracker to track the simulated targets using the radar detections. In this example, you learned how the tracker's history based logic promotes tracks. You also learned how the track uncertainty grows when a track is coasted and is reduced when the track is updated by a new detection.

Supporting Functions

initFilter

This function modifies the function initcvekf to handle higher velocity targets such as the airliners in the ATC 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 initialized track's 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 match filter's update rate
scaleAccelHorz = cast(1,classToUse);
scaleAccelVert = cast(1,classToUse);
Q = blkdiag(scaleAccelHorz^2, scaleAccelHorz^2, scaleAccelVert^2);
filter.ProcessNoise = Q;
end

tabulateData

This function returns a table comparing the ground truth and tracks

function truthTrackTable = tabulateData(scenario, tracks)
% Process truth data
platforms = scenario.Platforms(2:end); % Platform 1 is the radar
numPlats = numel(platforms);
trueAlt = zeros(numPlats,1);
trueSpd = zeros(numPlats,1);
trueHea = zeros(numPlats,1);
for i = 1:numPlats
    traj = platforms{i}.Trajectory;
    waypoints = waypointInfo(traj);
    times = waypoints.TimeOfArrival;
    waypoints = waypoints.Waypoints;
    trueAlt(i) = -waypoints(end,3);
    trueVel = (waypoints(end,:) - waypoints(end-1,:)) / (times(end)-times(end-1));
    trueSpd(i) = norm(trueVel) * 3600 / 1000; % Convert to km/h
    trueHea(i) = atan2d(trueVel(1),trueVel(2));
end
trueHea = mod(trueHea,360);

% Associate tracks with targets
atts = [tracks.ObjectAttributes];
tgtInds = [atts.TargetIndex];

% Process tracks assuming a constant velocity model
numTrks = numel(tracks);
estAlt = zeros(numTrks,1);
estSpd = zeros(numTrks,1);
estHea = zeros(numTrks,1);
truthTrack = zeros(numTrks,7);
for i = 1:numTrks
    estAlt(i) = -round(tracks(i).State(5));
    estSpd(i) = round(norm(tracks(i).State(2:2:6)) * 3600 / 1000); % Convert to km/h;
    estHea(i) = round(atan2d(tracks(i).State(2),tracks(i).State(4)));
    estHea(i) = mod(estHea(i),360);
    platID = tgtInds(i);
    platInd = platID - 1;
    truthTrack(i,:) = [tracks(i).TrackID, trueAlt(platInd), estAlt(i), trueHea(platInd), estHea(i), ...
        trueSpd(platInd), estSpd(i)];
end

% Organize the data in a table
names = {'TrackID','TrueAlt','EstimatedAlt','TrueHea','EstimatedHea','TrueSpd','EstimatedSpd'};
truthTrackTable = array2table(truthTrack,'VariableNames',names);
truthTrackTable = mergevars(truthTrackTable, (6:7), 'NewVariableName', 'Speed', 'MergeAsTable', true);
truthTrackTable.(6).Properties.VariableNames = {'True','Estimated'};
truthTrackTable = mergevars(truthTrackTable, (4:5), 'NewVariableName', 'Heading', 'MergeAsTable', true);
truthTrackTable.(4).Properties.VariableNames = {'True','Estimated'};
truthTrackTable = mergevars(truthTrackTable, (2:3), 'NewVariableName', 'Altitude', 'MergeAsTable', true);
truthTrackTable.(2).Properties.VariableNames = {'True','Estimated'};
truthTrackTable.TrackID = "T" + string(truthTrackTable.TrackID);
end