Create Scenario and Add Terrain

First create an Earth-centered tracking scenario, then add terrain data using the scenario object function groundSurface. You can specify the terrain as a matrix of height values, or like in this example, as a DTED file. The DTED used in this scenario spans latitude between 39 and 40 degrees North and longitude between 105 and 106 degrees West. This corresponds to a mountainous region in Colorado, USA.

scene = trackingScenario(IsEarthCentered=true, UpdateRate=1);
% Add terrain 
dtedfile = "n39_w106_3arc_v2.dt1";
groundSurface(scene,Terrain=dtedfile);

All surfaces added to the scenario are managed by the SurfaceManager object.

scene.SurfaceManager

ans = 
  SurfaceManager with properties:

    UseOcclusion: 1
        Surfaces: [1×1 fusion.scenario.GroundSurface]

scene.SurfaceManager.Surfaces(1)

ans = 
  GroundSurface with properties:

            Terrain: "n39_w106_3arc_v2.dt1"
    ReferenceHeight: 0
           Boundary: [2×2 double]

Next, add three platforms to the scenario. The first platform is a drone flying at a constant altitude of 20 meters above the ground. Use the following strategy to define a trajectory following the terrain:

% Add a platform on the ground
llacoords = [39.7894, -105.6284, 0;...
    39.8012, -105.6251, 0;...
    39.80449, -105.6420, 0];
groundTraj = geoTrajectory(llacoords,[0;120;240]);

numTrajectorySamples = 40;
toa = linspace(0,240, numTrajectorySamples);
pos = lookupPose(groundTraj, toa);
flightAltitude = 20; % 20 meters
pos(:,3) = height(scene.SurfaceManager,pos(:,[1 2])') + flightAltitude;

% Build the final trajectory and add platform
groundTraj = geoTrajectory(pos,toa);
drone = platform(scene,Trajectory=groundTraj);

The second platform is a ground vehicle that follows a mountain pass road. Coordinates and time values are saved in the roadCoordinates.mat file.

% Add a vehicle following the road
load roadCoordinates.mat
roadTraj = geoTrajectory(roadlla(10:end-30,:), 1.3*toas(10:end-30));
roadPlat = platform(scene,Trajectory=roadTraj);

The third platform is a radar tower. The tower is located on top of a mountain and the ground and aerial surveillance radar, mounted on it, stares at the ground and sky towards the Southeast. Set the HasElevation property to true to report elevation, which is important when tracking over terrain where elevation often varies.

% Add a tower
towerLatLon = [39.8057, -105.6246];
towerGroundHeight = height(scene.SurfaceManager,towerLatLon');
tower = platform(scene,Position=[ towerLatLon, towerGroundHeight]);
tower.Sensors = fusionRadarSensor(ScanMode="No scanning", ...
    UpdateRate=1/2,...
    SensorIndex=2,...
    FieldOfView=[180;35],...
    HasElevation=true,...
    ElevationResolution=1,...
    RangeResolution=5,...
    RangeLimits=[0 2000],...
    ReferenceRange=2000,...
    FalseAlarmRate=1e-7,...
    HasFalseAlarms=true,...
    MountingAngles=[200 -15 0],...
    MountingLocation=[0 0 -25],...
    HasINS=true, DetectionCoordinates='scenario',...
    HasNoise=true);

Define Tracking System

Use a multi-object tracker to track the drone and the ground vehicle. Use the default constant velocity motion model, which works reasonably well for tracking targets that move slowly. Increase slightly the AssignmentThreshold property to account for large measurement noise when the radar is detecting at long ranges.

towerTracker = trackerGNN(AssignmentThreshold=50);

Create Visualization and Simulate the Scenario

Visualize the scene using a trackingGlobeViewer object. By default, the viewer object does not display terrain. First, use the addCustomTerrain function to add the DTED file. Then specify the Terrain property as the DTED file name to visualize the terrain in the viewer. The addCustomTerrain function saves the DTED to a location and the function should be used only once. Also, specify the CoverageRangeScale property to scale down the radar coverage and reduce visual clutter. This, however, causes that the coverage displayed in the viewer no longer reflects the actual range of the radar. Record and visualize the occlusion history by using the occlusion object function of the SurfaceManager object. The legend used in the viewer is shown below.

try
    addCustomTerrain("TrackingOverTerrainExample",dtedfile);
catch
    disp("Custom Terrain was already added")
end
% Create and setup the trackingGlobeViewer object
f=uifigure;
globe = trackingGlobeViewer(f,Terrain="TrackingOverTerrainExample",CoverageRangeScale=0.1,TrackHistoryDepth=5e5);
campos(globe, [39.7861 -105.6287 3265]);
camorient(globe, [19.35 -23.18 0]);
plotPlatform(globe, tower, Marker='d');
plotCoverage(globe, coverageConfig(scene), "ECEF", Alpha=0.2);


% Initialize some variables for the simulation
towerRadarLLA = tower.Position + [0 0 25];
droneOcc = [];
groundOcc = [];
droneTracks = objectTrack.empty;
towerTracks = objectTrack.empty;

% Metric
ospa = trackOSPAMetric(Metric="OSPA",Distance="posabserr");
ospalog = [];

% Set random seed for repeatable results
s = rng(2022);

while advance(scene)
    time = scene.SimulationTime;

    % Generate detections
    [towerDets,~, config] = detect(tower, time);

    % Update tower tracker
    if config.IsValidTime && (isLocked(towerTracker) || ~isempty(towerDets))
        towerTracks = towerTracker(towerDets, time);
    elseif isLocked(towerTracker)
        towerTracks = predictTracksToTime(towerTracker, 'confirmed',time);
    end

    % Update globe
    plotPlatform(globe,[drone roadPlat],TrajectoryMode="Full",LineWidth=1);
    plotDetection(globe, towerDets, "ECEF");
    plotTrack(globe, towerTracks, "ECEF", LabelStyle="ATC",LineWidth=2);
    
    % Move camera
    moveGlobeCamera(globe, time);

    % Compute metric and save occlusion state for plots
    ospalog(end+1) = ospa(towerTracks,platformPoses(scene)); %#ok<SAGROW>
    droneOcc(end+1) = occlusion(scene.SurfaceManager,towerRadarLLA,drone.Position); %#ok<SAGROW> 
    groundOcc(end+1) = occlusion(scene.SurfaceManager,towerRadarLLA,roadPlat.Position); %#ok<SAGROW> 

    drawnow;
end

In the figure below, the occlusion plot on the left axis helps analyze the scenario by showing the occlusion status of the target over time. This provides the information to understand why radar detections are not available as well as why tracks are dropped during some periods. The OSPA metric plot on the right axis gives a quantitative assessment of the tracker performance and shows how tracking performance correlates with the occlusion status.

plotOcclusionAndMetric(droneOcc, groundOcc, ospalog);

The two GIFs below provide analysis of two specific periods in the simulation.

The first GIF above shows a segment of the simulation around 100 seconds, which corresponds to the time that the drone enters the occluded area at the foot of the radar-located mountain. Notice that the track is coasted, and its associated uncertainty grows due to the lack of observations until it is deleted. The previous OSPA plot slowly increases during coasting before jumping to the threshold value upon track deletion. In the meantime, the ground vehicle, travelling on the road, is occluded and therefore undetected.

The second GIF shows a later period of the simulation around 155 seconds. Both drone and ground vehicle are not occluded, as shown in the occlusion state plot, except for a moment when the drone passes by the saddle point in between mountains. The track is briefly coasted and recovered with the next available radar detections. This period of the simulation is also noticeable on the occlusion and OSPA plot. The performance starts to degrade (OSPA value increases) after the occlusion status switches but recovers immediately when the status switches again.

Last, take a snapshot of the globe viewer at the end of the simulation to get the full picture of the scenario.

campos(globe, [39.805195 -105.64668  3123]);
camorient(globe, [110 -5 0]);
drawnow;
figure;
snapshot(globe);

The true trajectories of the ground target and the drone are displayed in white whereas tracks are represented by colored lines. The three segments of the drone track are shown in yellow, green, and purple, respectively. The ground vehicle tracks are shown in blue and orange. In this particular scenario, the long occlusion time makes it difficult for the tracker to maintain a unique track ID for each target.

Clean Up

Reset random generator seed and remove the added DTED file.

rng(s);

% Remove custom terrain
if isvalid(f)
    close(f);
end
removeCustomTerrain("TrackingOverTerrainExample");

Conclusion

In this example you learned how to include terrain data from DTED files in a tracking scenario and how to use the SurfaceManager property of the trackingScenario to query height and occlusion information over the terrain. This allows to create trajectories that follow the terrain, such as the trajectory of a vehicle following a mountain pass road or a drone flying at constant altitude above ground. You also simulated radar detections that accounts for terrain occlusion and used a simple tracking system to track the targets. For targets that are difficult to recognize after a long time of terrain occlusion, alternate techniques relying on appearance or radar signal features could be used to re-identify lost vehicles.

Supporting functions

plotOcclusionAndMetric generates the plots of occlusion status between the tower and each platform as well as the track OSPA metric.

function plotOcclusionAndMetric(droneOcc, groundOcc, ospalog)
figure;
timevals = 1:numel(groundOcc);

yyaxis right
plot(timevals,ospalog, DisplayName="OSPA",LineWidth=3)
ylim([0 40]);
ylabel("OSPA");

yyaxis left
plot(timevals,droneOcc,DisplayName="Drone",LineWidth=3);
hold on
plot(timevals,groundOcc,Color=[0.2464    0.1569    0.6847],...
    DisplayName="Ground Target",...
    LineStyle='-',LineWidth=3);
yticklabels({'Visible','Occluded'});
yticks([0 1]);
ylim([-0.2 1.5]);

title("Occlusion Status and Tracker Performance");
xlabel("time")
legend();
drawnow
end

moveGlobeCamera moves the camera on the globe at predefined times.

function moveGlobeCamera(globe, time)
if time == 50
% Move camera down the pass
    campos(globe, [39.7998 -105.60964 3016]);
    camorient(globe, [270 1 0]);
  
elseif time == 90
% Move camera all the way down to observe occlusion
    campos(globe, [39.7963 -105.6188 2936]);
    camorient(globe, [306 6 0]);
elseif time == 120
% Move camera up the road and raise view angle
    campos(globe, [39.7953 -105.6323 3061]);
    camorient(globe, [14 -6 0]);

elseif time == 165
% Move camera to last section of the road
    campos(globe, [39.79876 -105.6428 3166]);
    camorient(globe, [35 -9 0]);
    
end
end