Simulate and Track En-Route Aircraft in Earth-Centered Scenarios
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 .
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:
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
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
