Simulate, Detect, and Track Anomalies in a Landing Approach

R2026b

This example shows how to automatically detect deviations and anomalies in aircraft making final approaches to an airport runway. You model an ideal landing approach trajectory, generate variants from it, simulate tracks, and issue warnings as soon as the tracks deviate from safe landing rules.

Introduction

Landing is a safety-critical stage of flight. An aircraft in final approach to landing must align itself with the runway, gradually descend to the ground, and reduce its ground speed while keeping it safely above stall speed. All these steps are done to ensure that the aircraft touches the ground softly to reduce the risk to passengers and to avoid physical damage to the aircraft or the runway. These rules can be easily defined by an aviation professional or they can be inferred from tracking data using machine learning [1]. In this example, you assume that the rules are already defined.

Major airports usually have multiple runways oriented in different directions. Approaching aircraft are guided by the air traffic controllers in the airport tower to land on one of the runways that is best aligned against the direction of wind at that time. During the approach, controllers monitor the aircraft using tracking systems. Airport traffic has increased over the last decades and, with it, the workload on air traffic controllers has increased as well. As a result, there is a need to automatically and reliably alert air traffic controllers to aircraft that approach the landing point in an unsafe manner: not aligned well with the runway, descending too quickly or too slowly, or approaching too quickly or too slowly.

Define Landing Rules

You define a landing approach trajectory into Logan International Airport in Boston, MA using a geoTrajectory object. The trajectory waypoints are aligned with runway 22L, which runs from north-east to south-west, and a glide slope of 3 degrees. The time of arrival and climb rate are defined to slow the approaching aircraft to a safe speed and a smooth touchdown. Note that the positive value of climb rate is used for a descending trajectory. You use the trackingGlobeViewer to visualize the trajectory on the map.

baselineApproachTrajectory = geoTrajectory([42.7069 -70.8395 1500; 42.5403 -70.9203 950; 42.3736 -71.001 0],[0;180;400],...
    'ClimbRate', [6; 3.75; 3.75]);
fig = uifigure('Units','normalized','Position',[0.1 0.1 0.8 0.8]);
viewer = trackingGlobeViewer(Parent=fig);
viewer.PlatformHistoryDepth = 0;
viewer.TrackHistoryDepth = 0;
viewer.TrackLabelScale = 2;
campos(viewer, [42.3072 -70.8463 12455]);
camorient(viewer, [335 -34 0]);
plotTrajectory(viewer, baselineApproachTrajectory, Color=[15 255 255]/255, LineWidth=1);

For a trajectory coming to land at Logan airport on runway 22L to be safe, the trajectory must satisfy the following rules:

• The trajectory must be closely aligned with the runway direction.

• The glide slope must be between 2.5 and 4 degrees in the last 20963 meters. At distances above 20963 meters, the altitude must be at least 3000 ft.

• The speed must be between 120 knots and 180 knots at the landing point. The upper speed bound can increase linearly with distance from the landing point.

You define these rules using the helperTrajectoryValidationRule (see the supporting file) and visualize the rules on the map using the helperShowRules.

% Define and show trajectory rules.
trajRules = defineTrajectoryRules();
helperShowRules(viewer, trajRules);
pause(2);
img = snapshot(viewer);
f = figure(Units="pixels",Position=[100 100 size(img,2) size(img,1)]);
ax = axes(f);
imshow(img,Parent=ax);

Define Scenario

In this section, you define an Earth-centered tracking scenario to simulate the trajectories of 20 incoming aircraft. Aircraft approaching landing in an airport are scheduled to avoid aerodynamic impact from one aircraft on the one following in its wake. The minimum safe time difference between two aircraft is one minute.

% Create an Earth-centered tracking scenario
scenario = trackingScenario('UpdateRate', 1, 'IsEarthCentered', true);

You use the perturbations object function to define a normal distribution around the baseline trajectory, baselineApproachTrajectory. Each waypoint in the trajectory is perturbed with a zero-mean normal distribution and a standard deviation that gets smaller from the first waypoint to the last (the landing point). At the first waypoint, the standard deviation is 5e-3 degrees in longitude and 300 meters in altitude. The standard deviation decreases to 1e-3 degrees in longitude and 150 meters in altitude in the mid-point and then to 1e-4 degrees in longitude and 0 in altitude at the end-point on the ground.

% Define perturbation to the approach trajectory
perturbations(baselineApproachTrajectory, 'Waypoints', 'Normal', zeros(3,3), [0 5e-3 300; 0 1e-3 150; 0 1e-4 0]);

To create 20 trajectories that are perturbed from the baseline trajectory, first clone the trajectory and then perturb it to apply perturbations to the trajectory waypoints. Then, you use the perturbations and perturb object functions again to perturb the TimeOfArrival of each trajectory and to make sure no additional perturbations are applied to the Waypoints. Then you attach each trajectory to a new platform. To perturb the entire scenario, you use the perturb object function.

% Generate perturbed trajectories
s = rng(2021, 'twister'); % Set random noise generator for repeatable results
numTrajectories = 20;
trajectories = cell(1,numTrajectories);
for i = 1:numTrajectories
    % Create a new trajectory with perturbed waypoints
    newTrajectory = clone(baselineApproachTrajectory);
    perturb(newTrajectory);
    
    % Perturb time of arrival of the new trajectory
    perturbations(newTrajectory, 'TimeOfArrival', 'Uniform',(i-1)*60, (i-1)*60+10);
    perturbations(newTrajectory, 'Waypoints', 'None');
    
    % Attach the trajectory to the platform
    platform(scenario, 'Trajectory', newTrajectory);

    % Store trajectories for analysis
    trajectories{i} = newTrajectory;
end

perturb(scenario);

To see which perturbed trajectories satisfy the rules of a safe approach to landing, you use the helper function validateTrajectory, provided at the bottom of this page. The function declares a trajectory to be anomalous if at least 1% of the points sampled from it violate any trajectory rule.

[truthAnomalyFlags, truthPercentAnomalousSteps] = validateTrajectory(trajectories, trajRules);

Plot the trajectories in yellow for anomalous trajectories and cyan for safe approaches. Overall, there are 7 anomalous trajectories out of the 20 generated trajectories.

plotTrajectory(viewer, trajectories(truthAnomalyFlags), Color=[255 255 17]/255, LineWidth=1);
plotTrajectory(viewer, trajectories(~truthAnomalyFlags), Color=[15 255 255]/255, LineWidth=1);
campos(viewer, [42.4808 -70.916 1136]);
camorient(viewer, [340 0 0]);
imshow(snapshot(viewer),Parent=ax);

Simulate Tracking System

Like other major airports in the USA, Logan uses an Airport Surface Detection Equipment - Model X (ASDE-X) to track aircraft during final approach and on the ground [2]. ASDE-X relies on an airport surveillance radar, automatic dependent surveillance-broadcast (ADS-B) reports from the approaching aircraft, and other methods to provide accurate tracking which is updated every second (for more details, see [1]).

To simplify the model of this tracking system, you model the airport tracking system using a statistical radar model, fusionRadarSensor, attached to the airport tower, and connect it to a Global Nearest Neighbor (GNN) tracker. Radar detections are usually imperfect; measurement noise corrupts the reported position of each target, the radar may generate false alarms due to clutter, and the radar may miss targets due to imperfect detection algorithms. These imperfections are addressed with the GNN tracker, which associates radar detections over time, and helps filter out clutter and maintain smoother trajectory estimates.

asdex = fusionRadarSensor(1, ...
    'ScanMode', 'No Scanning', ...
    'MountingAngles', [0 0 0], ...
    'FieldOfView', [360;20], ...
    'UpdateRate', 1, ...
    'ReferenceRange', 40000,...
    'RangeLimits', [0 50000], ...
    'RangeResolution', 100, ...
    'HasElevation', true, ...
    'HasRangeRate', false, ...
    'DetectionCoordinates', 'Sensor spherical', ...
    'HasINS',true,...
    'FalseAlarmRate', 1e-7, ...
    'ElevationResolution', 0.4, ...
    'AzimuthResolution', 0.4);
platform(scenario, 'Position', [42.3606 -71.011 0], 'Sensors', asdex);

Next, you define a multi-object GNN tracker using the task-oriented tracking approach provided in the toolbox. First, you define a target specification for a passenger aircraft modeled in geographic coordinates. Then, you define a sensor specification for a monostatic radar that matches the configuration of the simulated radar. For more details on how to configure and use the tracker, refer to the Track Heterogeneous Aircraft Using Heterogeneous Sensors.

% Define target specification
targetSpec = trackerTargetSpec('aerospace', 'aircraft', 'passenger');
targetSpec.IsGeographic = true;

% Define sensor specification matching the radar parameters
sensorSpec = trackerSensorSpec('aerospace', 'radar', 'monostatic');
sensorSpec.IsGeographic = true;
sensorSpec.PlatformPosition = [42.3606 -71.011 0];
sensorSpec.FieldOfView = [360 20];
sensorSpec.RangeLimits = [0 50000];
sensorSpec.RangeRateLimits = [-500 500];
sensorSpec.RangeResolution = 100;
sensorSpec.AzimuthResolution = 0.4;
sensorSpec.ElevationResolution = 0.4;
sensorSpec.FalseAlarmRate = 1e-7;
sensorSpec.DetectionProbability = 0.9;
sensorSpec.HasRangeRate = false;
sensorSpec.RangeRateResolution = 10;

Next, you create the tracker using the multiSensorTargetTracker with the GNN algorithm. You use a maximum Mahalanobis distance of 8 units to gate measurements further away from the track's prediction. You also confirm a track if it is reported by the radar in 4 out of 5 tracker updates.

tracker = multiSensorTargetTracker(targetSpec, sensorSpec, 'gnn');
tracker.MaxMahalanobisDistance = 8;
tracker.NumUpdatesForConfirmation = [4 5];

Then, you create assignment metrics to assign tracks to simulated truths. The assignment between tracks and truth is used to analyze the accuracy of the anomaly classification.

tam = trackAssignmentMetrics('AssignmentThreshold', 100, 'DivergenceThreshold', 200);

Run Scenario and Detect Anomalies

In this section, you run the full simulation and classify tracked objects to generate anomaly warnings. At each step, detections from the radar are converted to the tracker's expected data format and passed to the tracker. You use the validateTracks helper function to generate anomaly warnings for the tracks. Tracks that violate the safe approach rules are shown in yellow while tracks that follow these rules are shown in cyan.

% Clean the display and prepare it for the simulation.
clear(viewer)
campos(viewer, [42.3072 -70.8463 12455]);
camorient(viewer, [335 -34 0]);
helperShowRules(viewer, trajRules)
clear validateTracks
restart(scenario);

% Main loop
while advance(scenario)
    % Collect detections
    dets = detect(scenario);

    % Convert detections to tracker data format and update tracker
    sensorData = convertDetections(dets, scenario.SimulationTime);
    tracks = tracker(sensorData);
   
    % Get platform poses and assignment between tracks and truths.
    poses = platformPoses(scenario,"Quaternion","CoordinateSystem","Cartesian");
    tam(tracks,poses);
    [assignedTrackIDs, assignedTruthIDs] = currentAssignment(tam);

    % Validate the tracks with rules to find anomalous tracks.
    [tracksClassified, trackAnomalyHistory] = validateTracks(tracks, trajRules, assignedTrackIDs, assignedTruthIDs);

    % Visualize
    plotPlatform(viewer,scenario.Platforms,'ECEF',LabelStyle="None");
    plotSensorData(viewer,sensorData,sensorSpec,Color=[255 19 166]/255);
    colors = getTrackColorsByType(tracksClassified);
    plotTrack(viewer,tracks,targetSpec,Color=colors);

    % Take a snapshot at 930 seconds
    if abs(scenario.SimulationTime - 930) < 0.01
        imshow(snapshot(viewer),Parent=ax);
    end
end

The picture above shows the snapshot at 930 seconds into the simulation. Note that tracks identified as safe are colored cyan and tracks identified as anomalous are colored yellow. This identification is done at every simulation step as can be seen for track 1819, colored in yellow and for track 1993, colored in cyan.

Analyze Classification

To verify that the anomaly warnings were issued for the right tracks, you use the analyze helper function, shown at the bottom of this page.

The function uses the trackAnomalyHistory collected during the simulation and compares it to the truthPercentAnomalousSteps calculated for each trajectory. Similar to truth, tracks are assigned anomaly flags if they were declared anomalous at least 1% of the time steps. You can see that the anomalies are issued correctly for the seven trajectories that were found to be anomalous, with one false anomaly for TruthID = 9.

comparisonTable = analyze(trackAnomalyHistory ,truthPercentAnomalousSteps);
disp(comparisonTable)

    TruthID    Truth Anomaly Flag    Track Anomaly Flag
    _______    __________________    __________________

       1             false                 false       
       2             false                 false       
       3             false                 false       
       4             false                 false       
       5             false                 false       
       6             false                 false       
       7             true                  true        
       8             true                  true        
       9             false                 true        
      10             true                  true        
      11             false                 false       
      12             false                 false       
      13             true                  true        
      14             true                  true        
      15             false                 false       
      16             false                 false       
      17             false                 false       
      18             false                 false       
      19             true                  true        
      20             true                  true        

confusionchart(comparisonTable.("Truth Anomaly Flag"),comparisonTable.("Track Anomaly Flag"));

Summary

In this example, you learned how to use tracking data to generate real-time warnings for anomalies like an unsafe landing approach.

You used geoTrajectory to define an ideal landing approach trajectory in geographic coordinates. You then used perturbations and perturb to create 20 trajectories that deviate from the ideal landing approach trajectory and to schedule the trajectories one after the other in the trackingScenario. To model an airport tracking system, you used a fusionRadarSensor System object to simulate radar detections and a multiSensorTargetTracker to define a task-oriented GNN tracker configured to track passenger aircraft.

References

1. Raj Deshmukh and Inseok Hwang, "Anomaly Detection Using Temporal Logic Based Learning for Terminal Airspace Operations", AIAA SciTech Forum, 2019.

2. Federal Aviation Administration, "Fact Sheet - Airport Surface Detection Equipment, Model X (ASDE-X)". Retrieved May 2020.

Supporting Functions

convertDetections Convert objectDetection cell array to tracker data format

function sensorData = convertDetections(dets, time)
sensorData.LookTime = time;
sensorData.LookAzimuth = 0;
sensorData.LookElevation = 0;
numDets = numel(dets);
detTime = zeros(1, numDets);
az = zeros(1, numDets); el = az; rng = az;
azAcc = az; elAcc = az; rngAcc = az;
for d = 1:numDets
    detTime(d) = dets{d}.Time;
    az(d) = dets{d}.Measurement(1);
    el(d) = dets{d}.Measurement(2);
    rng(d) = dets{d}.Measurement(3);
    azAcc(d) = sqrt(dets{d}.MeasurementNoise(1,1));
    elAcc(d) = sqrt(dets{d}.MeasurementNoise(2,2));
    rngAcc(d) = sqrt(dets{d}.MeasurementNoise(3,3));
end
sensorData.DetectionTime = detTime;
sensorData.Azimuth = az;
sensorData.Elevation = el;
sensorData.Range = rng;
sensorData.AzimuthAccuracy = azAcc;
sensorData.ElevationAccuracy = elAcc;
sensorData.RangeAccuracy = rngAcc;
end

defineTrajectoryRules Define trajectory rules

function trajRules = defineTrajectoryRules
% This function defines rules for safe approach to landing on runway 22L at
% Logan International Airport in Boston, MA.

% The function uses the helperTrajectoryValidationRule attached as a
% supporting file to this example

% The trajectory must be closely aligned with the runway direction.
longitudeRule = helperTrajectoryValidationRule([42.37 42.71], [0.4587, -90.4379], [0.5128, -92.730]);

% The glide slope must be between 2.5 and 4 degrees in the last 20963
% meters. At distances above 20963 meters, the altitude must be at least
% 3000 ft. The rules are relative to range from the runway landing point.
altitudeRule1 = helperTrajectoryValidationRule([100 20963], [sind(2.5) 0], [sind(4) 0]);
altitudeRule2 = helperTrajectoryValidationRule([20963 40000], 3000 * 0.3048, [sind(4) 0]);

% The speed must be between 120 knots and 180 knots at the landing point.
% The upper speed bound can increase linearly with distance from the
% landing point.
speedRule = helperTrajectoryValidationRule([0 40000], 61.733, [1e-3 100]);

% Collect all the rules.
trajRules = [longitudeRule;altitudeRule1;altitudeRule2;speedRule];
end

validateTrajectory Validate each trajectory

function [truthAnomalyFlags, percentAnomalousSteps] = validateTrajectory(trajectories, rules)
numTrajectories = numel(trajectories);
numAnomalousSteps = zeros(1, numTrajectories);
for tr = 1:numTrajectories
    if iscell(trajectories)
        traj = trajectories{tr};
    elseif numTrajectories == 1
        traj = trajectories;
    else
        traj = trajectories(tr);
    end
    timesamples = (traj.TimeOfArrival(1):traj.TimeOfArrival(end));
    [pos,~,vel] = lookupPose(traj, timesamples);
    posECEF = lookupPose(traj, timesamples, 'ECEF');
    landingPoint = [1.536321  -4.462053   4.276352]*1e6;

    for i = 1:numel(timesamples)
        distance = norm(posECEF(i,:) - landingPoint);
        isLongitudeValid = validate(rules(1),pos(i,1),pos(i,2));
        isAltitudeValid = (validate(rules(2),distance,pos(i,3)) && validate(rules(3),distance,pos(i,3)));
        isSpeedValid = validate(rules(4),distance,norm(vel(i,:)));
        isValid = isLongitudeValid && isAltitudeValid && isSpeedValid;
        numAnomalousSteps(tr) = numAnomalousSteps(tr) + ~isValid;
    end
end
percentAnomalousSteps = numAnomalousSteps ./ numel(timesamples) * 100;
truthAnomalyFlags = (percentAnomalousSteps > 1);
end

validateTracks Validate tracks vs anomaly rules

function [tracks, history] = validateTracks(tracks, rules, assignedTrackIDs, assignedTruthIDs)
persistent trackAnomalyHistory
if isempty(trackAnomalyHistory)
    trackAnomalyHistory = repmat(struct('TrackID', 0, 'AssignedTruthID', 0, 'NumSteps', 0, 'NumAnomalousSteps', 0),30,1);
end

posECEF = getTrackPositions(tracks, [1 0 0 0 0 0; 0 0 1 0 0 0; 0 0 0 0  1 0]);
pos = fusion.internal.frames.ecef2lla(posECEF);
vel = getTrackVelocities(tracks, [0 1 0 0 0 0; 0 0 0 1 0 0; 0 0 0 0 0 1]);
numTracks = numel(tracks);
landingPoint = [1.536321  -4.462053   4.276352]*1e6;

trackIDs = [trackAnomalyHistory.TrackID];

for tr = 1:numTracks
    % Only validate tracks if altitude is greater than 0
    if pos(tr,3) > 15
        distance = norm(posECEF(tr,:) - landingPoint);
        isLongitudeValid = validate(rules(1),pos(tr,1),pos(tr,2));
        isAltitudeValid = (validate(rules(2),distance,pos(tr,3)) && validate(rules(3),distance,pos(tr,3)));
        isSpeedValid = validate(rules(4),distance,norm(vel(tr,:)));
        isValid = isLongitudeValid && isAltitudeValid && isSpeedValid;
    else
        isValid = true;
    end

    tracks(tr).ObjectClassID = uint8(~isValid) + uint8(isValid)*6; % To get the right color for tracks

    % Update anomaly history
    inHistory = (tracks(tr).TrackID == trackIDs);
    if any(inHistory)
        trackAnomalyHistory(inHistory).NumSteps = trackAnomalyHistory(inHistory).NumSteps + 1;
        trackAnomalyHistory(inHistory).NumAnomalousSteps = trackAnomalyHistory(inHistory).NumAnomalousSteps + ~isValid;
    else
        ind = find(trackIDs == 0, 1, 'first');
        trackAnomalyHistory(ind).AssignedTruthID = assignedTruthIDs(tracks(tr).TrackID == assignedTrackIDs);
        trackAnomalyHistory(ind).TrackID = tracks(tr).TrackID;
        trackAnomalyHistory(ind).NumSteps = 1;
        trackAnomalyHistory(ind).NumAnomalousSteps = ~isValid;
    end
end
history = trackAnomalyHistory;
end

analyze - Analyze the track anomaly history and compare it to the truth anomaly percentage

function comparisonTable = analyze(trackAnomalyHistory, percentTruthAnomalous)
trackAnomalyHistory = trackAnomalyHistory([trackAnomalyHistory.TrackID] > 0);
anomalousSteps = [trackAnomalyHistory.NumAnomalousSteps];
numSteps = [trackAnomalyHistory.NumSteps];
trackAssignedTruths = [trackAnomalyHistory.AssignedTruthID];
assignedTruths = unique(trackAssignedTruths);
numTrackAnomalousSteps = zeros(numel(assignedTruths),1);
numTrackSteps = zeros(numel(assignedTruths),1);

for i = 1:numel(assignedTruths)
    inds = (assignedTruths(i) == trackAssignedTruths);
    numTrackAnomalousSteps(i) = sum(anomalousSteps(inds));
    numTrackSteps(i) = sum(numSteps(inds));
end
percentTrackAnomalous = numTrackAnomalousSteps ./ numTrackSteps * 100;
trueAnomaly = (percentTruthAnomalous > 1)';
anomalyFlags = (percentTrackAnomalous > 1);
comparisonTable = table((1:20)',trueAnomaly, anomalyFlags,...
    'VariableNames', {'TruthID', 'Truth Anomaly Flag', 'Track Anomaly Flag'});
end

getTrackColorsByType Returns the track color indicating if the track is nominal or an anomaly.

function colors = getTrackColorsByType(tracks)
colorsDark = [255 255 17;...
    19 159 255;...
    255 105 41;...
    100 212 19;...
    183 70 255;...
    15 255 255;...
    255 19 166]/255;
trackColorInds = [tracks.ObjectClassID];
colors = colorsDark(trackColorInds,:);
end