Detect and Track LEO Satellite Constellation with Ground Radars

Since R2021a

This example uses:

This example shows how to import a Two-Line Element (TLE) file of a satellite constellation, simulate radar detections of the constellation, and track the constellation.

The task of populating and maintaining a catalog of space objects orbiting Earth is crucial in space surveillance. This task consists of several processes: detecting and identifying new objects and adding them to the catalog, updating known objects orbits in the catalog, tracking orbit changes throughout their lifetime, and predicting reentries in the atmosphere. In this example, we are study how to detect and track new satellites and add them to a catalog.

To guarantee safe operations in space and prevent collisions with other satellites or known debris, it important to correctly detect and track newly launched satellites. Space agencies typically share prelaunch information, which can be used to select a search strategy. A Low Earth Orbit (LEO) satellite search strategy consisting of fence-type radar systems is commonly used. A fence-type radar system searches a finite volume in space and detects satellites as they pass through its field of view. This strategy can detect and track newly launched constellation quickly. [1]

Importing a satellite constellation from a TLE file

Two-Line Element sets are a common data format to save orbital information of satellites. You can use the satelliteScenario object to import satellite orbits defined in a TLE file. By default, the imported satellite orbits are propagated using the SGP4 orbit propagation algorithm which provides good accuracy for LEO objects. In this example, these orbits provide with the ground truth to test the radar tracking system capability to detect newly launched satellites.

% Create a satellite scenario
satscene = satelliteScenario;
% Add satellites from TLE file.
tleFile = "leoSatelliteConstellation.tle";
constellation = satellite(satscene, tleFile);

Use the satellite scenario viewer to visualize the constellation.

play(satscene);

Simulating synthetic detections and track constellation

Modeling space surveillance radars

Define two stations with fan-shaped radar beams looking into space. The fans cut through the satellite orbits to maximize the number of detections. The radar stations located on North America form an East-West fence.

% First station coordinates in LLA
station1 = [48 -80 0];

% Second station coordinates in LLA
station2 = [50 -117 0];

Each station is equipped with a radar, which is modeled by using a fusionRadarSensor object. In order to detect satellites in the LEO range, the radar has the following requirements:

Detecting a 10 dBsm object up to 2000 km away
Resolving objects horizontally and vertically with a precision of 100 m at 2000 km range
Having a field of view of 120 degrees in azimuth and 30 degrees in elevation
Looking up into space

% Create fan-shaped monostatic radars
fov = [120;40];
radar1 = fusionRadarSensor(1,...
    'UpdateRate',0.1,... 10 sec
    'ScanMode','No scanning',...
    'MountingAngles',[0 90 0],... look up
    'FieldOfView',fov,... degrees
    'ReferenceRange',2000e3,... m
    'RangeLimits',  [0 2000e3], ... m
    'ReferenceRCS', 10,... dBsm
    'HasFalseAlarms',false,...
    'HasNoise', true,...
    'HasElevation',true,...
    'AzimuthResolution',0.03,... degrees
    'ElevationResolution',0.03,... degrees
    'RangeResolution',2000, ... m % accuracy ~= 2000 * 0.05 (m)
    'DetectionCoordinates','Sensor Spherical',...
    'TargetReportFormat','Detections');

radar2 = clone(radar1);
radar2.SensorIndex = 2;

Radar Processing Chain

In this example, several coordinate conversions and axes transformation are performed to properly run the radar tracking chain. The diagram below illustrates how the inputs defined above are transformed and passed to the radar.

In the first step, you calculate each satellite pose in the local radar station NED axes. You achieve this by first obtaining the ground station ECEF pose and converting the satellite position and velocity to the ECEF coordinates. The radar input is obtained by taking the differences between the satellite pose and the ground station pose and rotating the differences into ground station local NED axes. See the assembleRadarInputs supporting function for the implementation details.

In the second step, you add the required information to the detection object so that the tracker can operate with an ECEF state. You use the MeasurementParameters property in each object detection to achieve that purpose, as shown in the addMeasurementParams supporting function.

Defining a tracker

The radar models you defined above output detections. To estimate the satellite orbits, you use a tracker. The Sensor Fusion and Tracking Toolbox™ provides a variety of multi-object trackers. In this example, you choose a Joint Probabilistic Data Association (JPDA) tracker because it offers a good balance of tracking performance and computational cost.

You need to define a tracking filter for the tracker. You can use a lower fidelity model than SGP4, such as a Keplerian integration of the equation of motion, to track the satellite. Often, the lack of fidelity in the motion model of targets is compensated by measurement updates and incorporating process noise in the filter. The supporting function initKeplerUKF defines the tracking filter.

% Define the tracker
tracker = trackerJPDA('FilterInitializationFcn',@initKeplerUKF,...
    'HasDetectableTrackIDsInput',true,...
    'ClutterDensity',1e-40,...
    'AssignmentThreshold',1e4,...
    'DeletionThreshold',[5 8],...
    'ConfirmationThreshold',[5 8]);

Running the simulation

In the remainder of this example, you step through the scenario to simulate radar detections and track satellites. This section uses the trackingGlobeViewer for visualization. You use this class to display sensor and tracking data with uncertainty ellipses and show the true position of each satellite.

viewer = trackingGlobeViewer('ShowDroppedTracks',false, 'PlatformHistoryDepth',700);

% Define coverage configuration of each radar and visualize it on the globe
ned1 = dcmecef2ned(station1(1), station1(2));
ned2 = dcmecef2ned(station2(1), station2(2));
covcon(1) = coverageConfig(radar1,lla2ecef(station1),quaternion(ned1,'rotmat','frame'));
covcon(2) = coverageConfig(radar2,lla2ecef(station2),quaternion(ned2, 'rotmat','frame'));
plotCoverage(viewer, covcon, 'ECEF');

You first generate the entire history of the states of the constellation over 5 hours. Then, you simulate radar detections and generate tracks in a loop.

satscene.StopTime = satscene.StartTime + hours(5);
satscene.SampleTime = 10;
numSteps = ceil(seconds(satscene.StopTime - satscene.StartTime)/satscene.SampleTime);

% Get constellation positions and velocity over the course of the simulation
plats = repmat(...
    struct('PlatformID',0,'Position',[0 0 0], 'Velocity', [0 0 0]),...
    numSteps, 40);
for i=1:numel(constellation)
    [pos, vel] = states(constellation(i),'CoordinateFrame',"ECEF");
    for j=1:numSteps
        plats(j,i).Position = pos(:,j)';
        plats(j,i).Velocity = vel(:,j)';
        plats(j,i).PlatformID = i;
    end
end

% Initialize tracks and tracks log
confTracks = objectTrack.empty(0,1);
trackLog = cell(1,numSteps);

% Initialize radar plots
radarplt = helperRadarPlot(fov);

% Set random seed for reproducible results
s = rng;
rng(2020);
step = 0;
while step < numSteps
    time = step*satscene.SampleTime;
    step = step + 1;
    
    % Generate detections of the constellation following the radar
    % processing chain
    targets1 = assembleRadarInputs(station1, plats(step,:));
    [dets1,numdets1] = radar1(targets1, time);
    dets1 = addMeasurementParams(dets1,numdets1,station1);
    
    targets2 = assembleRadarInputs(station2, plats(step,:));
    [dets2, numdets2] = radar2(targets2, time);
    dets2 = addMeasurementParams(dets2, numdets2, station2);
    
    detections = [dets1; dets2];
    updateRadarPlots(radarplt,targets1, targets2 ,dets1, dets2)
    
    % Generate and update tracks
    detectableInput = isDetectable(tracker,time, covcon);
    if ~isempty(detections) || isLocked(tracker)
        [confTracks,~,~,info] = tracker(detections,time,detectableInput);
    end
    trackLog{step} = confTracks;

    % Update viewer
    plotPlatform(viewer, plats(step,:),'ECEF', 'Color',[1 0 0],'LineWidth',1);
    plotDetection(viewer, detections,'ECEF');
    plotTrack(viewer, confTracks, 'ECEF','Color',[0 1 0], 'LineWidth',3);
            
end

{"String":"Figure contains 2 axes objects. Axes object 1 with title Radar 1 Field of View contains 772 objects of type line. Axes object 2 with title Radar 2 Field of View contains 1052 objects of type line.","Tex":["Radar 1 Field of View","Radar 2 Field of View"],"LaTex":[]}

The plot above shows the orbits (blue dots) and the detections (red circles) from the point of view of each radar.

% Restore previous random seed state
rng(s);

figure;
snapshot(viewer);

Figure contains an axes object. The axes object contains an object of type image.

After 5 hours of tracking, about half the constellation is tracked successfully. Maintaining tracks with a partial orbit coverage is challenging since satellites can often stay undetected for long periods of time in this configuration. In this example, there are only two radar stations. Additional tracking stations are expected to generate better tracking performance. The assignment metrics, which evaluate the tracking performance by comparing between the true objects and the tracks, are shown below.

% Show Assignment metrics
truthIdFcn = @(x)[x.PlatformID];

tam = trackAssignmentMetrics('DistanceFunctionFormat','custom',...
    'AssignmentDistanceFcn',@distanceFcn,...
    'DivergenceDistanceFcn',@distanceFcn,...
    'TruthIdentifierFcn',truthIdFcn,...
    'AssignmentThreshold',1000,...
    'DivergenceThreshold',2000);

for i=1:numSteps
    % Extract the tracker and ground truth at the i-th tracker update
    tracks = trackLog{i};
    truths = plats(i,:);
    % Extract summary of assignment metrics against tracks and truths
    [trackAM,truthAM] = tam(tracks, truths);
end

% Show cumulative metrics for each individual recorded truth object
results = truthMetricsTable(tam);
results(:,{'TruthID','AssociatedTrackID','BreakLength','EstablishmentLength'})

ans=40×4 table
    TruthID    AssociatedTrackID    BreakLength    EstablishmentLength
    _______    _________________    ___________    ___________________

       1               52                5                 232        
       2               24                0                 594        
       3                4                0                  56        
       4               55               11                 492        
       5               18                0                 437        
       6               48                8                 811        
       7               21                0                 513        
       8               27                0                 661        
       9               39                0                1221        
      10               50                0                1504        
      11               43                0                1339        
      12               37                0                1056        
      13              NaN                0                1800        
      14              NaN                0                1800        
      15              NaN                0                1800        
      16              NaN                0                1800        
      ⋮

The table above lists 40 satellites in the launched constellation and shows the tracked satellites with associated track IDs. A track ID of value NaN indicates that the satellite is not tracked by the end of the simulation. This either means that the orbit of the satellite did not pass through the field of view of one of the two radars or the track of the satellite has been dropped. The tracker can drop a track due to an insufficient number of initial detections, which leads to a large uncertainty on the estimate. Alternately, the tracker can drop the track if the satellite is not re-detected soon enough, such that the lack of updates leads to divergence and eventually deletion.

Summary

In this example, you have learned how to use the satelliteScenario object from the Aerospace Toolbox™ to import orbit information from TLE files. You propagated the satellite trajectories using SGP4 and visualized the scenario using the Satellite Scenario Viewer. You learned how to use the radar and tracker models from the Sensor Fusion and Tracking Toolbox™ to model a space surveillance radar tracking system. The constructed tracking system can predict the estimated orbit of each satellite using a low fidelity model.

Supporting functions

initKeplerUKF initializes an Unscented Kalman filter using the Keplerian motion model. Set Alpha = 1, Beta = 0, and Kappa = 0 to ensure robustness of the unscented filter over long prediction period.

function filter = initKeplerUKF(detection)

radarsphmeas = detection.Measurement;
[x, y, z] = sph2cart(deg2rad(radarsphmeas(1)),deg2rad(radarsphmeas(2)),radarsphmeas(3));
radarcartmeas = [x y z];
Recef2radar = detection.MeasurementParameters.Orientation;
ecefmeas = detection.MeasurementParameters.OriginPosition + radarcartmeas*Recef2radar;
% This is equivalent to:
% Ry90 = [0 0 -1 ; 0 1 0; 1 0 0]; % frame rotation of 90 deg around y axis
% nedmeas(:) = Ry90' * radarcartmeas(:);
% ecefmeas = lla2ecef(lla) +  nedmeas * dcmecef2ned(lla(1),lla(2));
initState = [ecefmeas(1); 0; ecefmeas(2); 0; ecefmeas(3); 0];

sigpos = 2;% m
sigvel = 0.5;% m/s^2

filter = trackingUKF(@keplerorbit,@cvmeas,initState,...
    'Alpha', 1, 'Beta', 0, 'Kappa', 0, ...
    'StateCovariance', diag(repmat([1000, 10000].^2,1,3)),...
    'ProcessNoise',diag(repmat([sigpos, sigvel].^2,1,3)));

end

function state = keplerorbit(state,dt)
% keplerorbit performs numerical integration to predict the state of
% Keplerian bodies. The state is [x;vx;y;vy;z;vz]

% Runge-Kutta 4th order integration method:
k1 = kepler(state);
k2 = kepler(state + dt*k1/2);
k3 = kepler(state + dt*k2/2);
k4 = kepler(state + dt*k3);

state = state + dt*(k1+2*k2+2*k3+k4)/6;

    function dstate=kepler(state)
        x =state(1,:);
        vx = state(2,:);
        y=state(3,:);
        vy = state(4,:);
        z=state(5,:);
        vz = state(6,:);

        mu = 398600.4405*1e9; % m^3 s^-2
        omega = 7.292115e-5; % rad/s
        
        r = norm([x y z]);
        g = mu/r^2;
        
        % Coordinates are in a non-intertial frame, account for Coriolis
        % and centripetal acceleration
        ax = -g*x/r + 2*omega*vy + omega^2*x;
        ay = -g*y/r - 2*omega*vx + omega^2*y;
        az = -g*z/r;
        dstate = [vx;ax;vy;ay;vz;az];
    end
end

isDetectable is used in the example to determine which tracks are detectable at a given time.

function detectInput = isDetectable(tracker,time,covcon)

if ~isLocked(tracker)
    detectInput = zeros(0,1,'uint32');
    return
end
tracks = tracker.predictTracksToTime('all',time);
if isempty(tracks)
    detectInput = zeros(0,1,'uint32');
else
    alltrackid = [tracks.TrackID];
    isDetectable = zeros(numel(tracks),numel(covcon),'logical');
    for i = 1:numel(tracks)
        track = tracks(i);
        pos_scene = track.State([1 3 5]);
        for j=1:numel(covcon)
            config = covcon(j);
            % rotate position to sensor frame:
            d_scene = pos_scene(:) - config.Position(:);
            scene2sens = rotmat(config.Orientation,'frame');
            d_sens = scene2sens*d_scene(:);
            [az,el] = cart2sph(d_sens(1),d_sens(2),d_sens(3));
            if abs(rad2deg(az)) <= config.FieldOfView(1)/2 && abs(rad2deg(el)) < config.FieldOfView(2)/2
                isDetectable(i,j) = true;
            else
                isDetectable(i,j) = false;
            end
        end
    end
    
    detectInput = alltrackid(any(isDetectable,2))';
end
end

assembleRadarInput is used to construct the constellation poses in each radar body frame. This is the first step described in the diagram.

function targets = assembleRadarInputs(station, constellation)
% For each satellite in the constellation, derive its pose with respect to
% the radar frame.
% inputs:
%          - station        :  LLA vector of the radar ground station
%          - constellation  :  Array of structs containing the ECEF position
%                              and ECEF velocity of each satellite
% outputs:
%          - targets        :  Array of structs containing the pose of each
%                              satellite with respect to the radar, expressed in the local
%                              ground radar frame (NED)

% Template structure for the outputs which contains all the field required
% by the radar step method
targetTemplate =  struct( ...
                'PlatformID', 0, ...
                'ClassID', 0, ...
                'Position', zeros(1,3), ...
                'Velocity', zeros(1,3), ...
                'Acceleration', zeros(1,3), ...
                'Orientation', quaternion(1,0,0,0), ...
                'AngularVelocity', zeros(1,3),...
                'Dimensions', struct( ...
                             'Length', 0, ...
                             'Width', 0, ...
                             'Height', 0, ...
                             'OriginOffset', [0 0 0]), ...
                'Signatures', {{rcsSignature}});


% First fill in the current satellite ECEF pose
targetPoses = repmat(targetTemplate, 1, numel(constellation));
for i=1:numel(constellation)
    targetPoses(i).Position = constellation(i).Position;
    targetPoses(i).Velocity = constellation(i).Velocity;
    targetPoses(i).PlatformID = constellation(i).PlatformID;
    % Orientation and angular velocity are left null, assuming satellite to
    % be point targets with a uniform rcs
end

% Then derive the radar pose in ECEF based on the ground station location
Recef2station = dcmecef2ned(station(1), station(2));
radarPose.Orientation = quaternion(Recef2station,'rotmat','frame');
radarPose.Position = lla2ecef(station);
radarPose.Velocity = zeros(1,3);
radarPose.AngularVelocity = zeros(1,3);

% Finally, take the difference and rotate each vector to the ground station
% NED axes
targets = targetPoses;
for i=1: numel(targetPoses)
    thisTgt = targetPoses(i);
    pos = Recef2station*(thisTgt.Position(:) - radarPose.Position(:));
    vel = Recef2station*(thisTgt.Velocity(:) - radarPose.Velocity(:)) - cross(radarPose.AngularVelocity(:),pos(:));
    angVel = thisTgt.AngularVelocity(:) - radarPose.AngularVelocity(:);
    orient = radarPose.Orientation' * thisTgt.Orientation;

    % Store into target structure array
    targets(i).Position(:) = pos;
    targets(i).Velocity(:) = vel;
    targets(i).AngularVelocity(:) = angVel;
    targets(i).Orientation = orient;
end
end

addMeasurementParam implements step 2 in the radar chain process as described in the diagram.

function dets = addMeasurementParams(dets, numdets, station)
% Add radar station information to the measurement parameters
Recef2station = dcmecef2ned(station(1), station(2));
for i=1:numdets
    dets{i}.MeasurementParameters.OriginPosition = lla2ecef(station);
    dets{i}.MeasurementParameters.IsParentToChild = true; % parent = ecef, child = radar
    dets{i}.MeasurementParameters.Orientation = dets{i}.MeasurementParameters.Orientation' * Recef2station;
end
end

distanceFcn is used with the assignment metrics to evaluate the tracking assignment.

function d = distanceFcn(track, truth)

estimate = track.State([1 3 5 2 4 6]);
true = [truth.Position(:) ; truth.Velocity(:)];
cov = track.StateCovariance([1 3 5 2 4 6], [1 3 5 2 4 6]);
d = (estimate - true)' / cov * (estimate - true);
end

Reference

[1] Sridharan, Ramaswamy, and Antonio F. Pensa, eds. Perspectives in Space Surveillance. MIT Press, 2017.