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