Обнаружьте и отследите спутниковое созвездие LEO с наземными радарами

В этом примере показано, как импортировать файл 2D линейного элемента (TLE) спутникового созвездия, симулируйте радарные обнаружения созвездия и отследите созвездие.

Задача заполнения и поддержания каталога объектов пробела, вращающихся вокруг Земли, крайне важна для наблюдения за космическим пространством. Эта задача состоит из нескольких процессов: обнаружение и идентификация новых объектов и добавление их к каталогу, обновление известных орбит объектов в каталоге, отслеживание изменений орбиты в течение их времени жизни и предсказания возвращений в атмосфере. В этом примере мы - исследование, как обнаружить и отследить новые спутники и добавить их в каталог.

Гарантировать безопасную работу на пробеле и предотвратить столкновения с другими спутниками или известными обломками, это важный, чтобы правильно обнаружить и отследить недавно запущенные спутники. Космические агентства обычно делятся информацией перед запуском, которая может использоваться, чтобы выбрать поисковую стратегию. Стратегия поиска спутника Низкой околоземной орбиты (LEO), состоящая из радиолокационных систем типа забора, обычно используется. Радиолокационная система типа забора ищет конечный объем на пробеле и обнаруживает спутники, когда они проходят через его поле зрения. Эта стратегия может обнаружить и отследить недавно запущенное созвездие быстро. [1]

Импорт спутникового созвездия из файла TLE

Наборы 2D линейного элемента являются форматом общих данных, чтобы сохранить орбитальную информацию спутников. Можно использовать satelliteScenario возразите, чтобы импортировать спутниковые орбиты, заданные в файле TLE. По умолчанию импортированные спутниковые орбиты распространены с помощью алгоритма распространения орбиты SGP4, который обеспечивает хорошую точность для объектов LEO. В этом примере эти орбиты предоставляют основную истину, чтобы протестировать радарную возможность системы слежения обнаружить недавно запущенные спутники.

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

Используйте спутниковое средство просмотра сценария, чтобы визуализировать созвездие.

play(satscene);

Симуляция синтетических обнаружений и созвездия дорожки

Моделирование радаров наблюдения за космическим пространством

Задайте две станции с веерообразными радарными лучами, смотрящими в космос. Вентиляторы прорубают спутниковые орбиты, чтобы максимизировать количество обнаружений. Радарные станции, расположенные на Северной Америке, формируют забор восток - запад.

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

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

Каждая станция оборудована радаром, который моделируется при помощи fusionRadarSensor объект. Для того, чтобы обнаружить спутники в области значений LEO, радар имеет следующие требования:

  • Обнаружение 10 объектов dBsm на расстоянии в 2 000 км

  • Решение объектов горизонтально и вертикально с точностью 100 м в 2 000-километровой области значений

  • Наличие поля зрения 120 градусов в области азимута и 30 градусов в области вертикального изменения

  • Поиск в космос

% 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;

Радарная цепь обработки

В этом примере несколько координатных преобразований преобразований и осей выполняются, чтобы правильно запустить радарную цепь отслеживания. Схема ниже иллюстрирует, как входные параметры, заданные выше, преобразовываются и передаются радару.

В первом шаге вы вычисляете каждое спутниковое положение в локальной радарной станции оси NED. Вы достигаете этого первым получением наземной станции положение ECEF и преобразование спутникового положения и скорости к координатам ECEF. Радарный вход получен путем взятия различий между спутниковым положением и положением наземной станции и вращения различий в наземную станцию локальные оси NED. Смотрите assembleRadarInputs поддерживание функции для деталей реализации.

На втором шаге вы добавляете запрошенную информацию в объект обнаружения так, чтобы средство отслеживания могло действовать с состоянием ECEF. Вы используете MeasurementParameters свойство в каждом обнаружении объектов, чтобы достигнуть той цели, как показано в addMeasurementParams поддерживание функции.

Определение средства отслеживания

Радарные модели вы задали выше выходных обнаружений. Чтобы оценить спутниковые орбиты, вы используете средство отслеживания. Sensor Fusion and Tracking Toolbox™ обеспечивает множество мультиобъектных средств отслеживания. В этом примере вы выбираете средство отслеживания Joint Probabilistic Data Association (JPDA), потому что он предлагает хороший баланс отслеживания эффективности и вычислительной стоимости.

Необходимо задать фильтр отслеживания для средства отслеживания. Можно использовать более низкую модель точности, чем SGP4, такой как Кеплеровское интегрирование уравнения движения, чтобы отследить спутник. Часто, отсутствие точности в модели движения целей компенсируется обновлениями измерения и соединяющимся шумом процесса в фильтре. Функция поддержки initKeplerUKF задает фильтр отслеживания.

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

Выполнение симуляции

В оставшейся части этого примера вы продвигаетесь через сценарий, чтобы симулировать радарные обнаружения и спутники дорожки. Этот раздел использует отдельный класс помощника helperScenarioGlobeViewer для визуализации. Можно использовать этот класс помощника, чтобы отобразить датчик и данные об отслеживании с замещающими знаками неопределенности и показать истинное положение каждого спутника.

globeDisplay = helperScenarioGlobeViewer;

% 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(globeDisplay, covcon)

Вы сначала генерируете целую историю состояний созвездия более чем 5 часов. Затем вы симулируете радарные обнаружения и генерируете дорожки в цикле.

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',"Geographic");
    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);
%         confTracks = deleteBadTracks(tracker,confTracks);
    end
    trackLog{step} = confTracks;

    % Update globe display
    updateDisplay(globeDisplay,time,plats(step,:),detections,[],confTracks);
            
end

График выше показов орбиты (синие точки) и обнаружения (красные круги) с точки зрения каждого радара.

% Restore previous random seed state
rng(s);
figure;
snap(globeDisplay);

После 5 часов отслеживания приблизительно половина созвездия прослежена успешно. Поддержание дорожек с частичным покрытием орбиты сложно, поскольку спутники могут часто оставаться необнаруженными в течение долгих промежутков времени в этой настройке. В этом примере существует только две радарных станции. Дополнительные станции отслеживания, как ожидают, сгенерируют лучше отслеживающую эффективность. Метрики присвоения, которые оценивают эффективность отслеживания путем сравнения между истинными объектами и дорожками, показывают ниже.

% 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,:);
    % For correct truth-track assignment evaluation, the truth is converted
    % to ECEF
    for j=1:numel(truths)
        R = dcmecef2ned(truths(j).Position(1), truths(j).Position(2));
        truths(j).Velocity = truths(j).Velocity * R;
        truths(j).Position = lla2ecef(truths(j).Position);
    end
    % 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               57                5                 232        
       2               46                0                1283        
       3               32                0                 668        
       4               62               12                 492        
       5               22                0                 436        
       6               54                5                 807        
       7               26                0                 513        
       8               31                0                 665        
       9               43                0                1221        
      10               56                0                1500        
      11               49                0                1339        
      12               41                0                1056        
      13              NaN                0                1800        
      14              NaN                0                1800        
      15              NaN                0                1800        
      16              NaN                0                1800        
      ⋮

Таблица выше приводит 40 спутников в запущенном созвездии и показывает отслеженные спутники со связанными идентификаторами дорожки. ID дорожки значения, NaN указывает, что спутник не прослежен к концу симуляции. Это или означает, что орбита спутника не проходила через поле зрения одного из этих двух радаров или дорожки спутника, был пропущен. Средство отслеживания может пропустить дорожку из-за недостаточного количества начальных обнаружений, которое приводит к большой неопределенности на оценке. Альтернативно, средство отслеживания может пропустить дорожку, если спутник не повторно обнаруживается достаточно скоро, такой, что отсутствие обновлений приводит к расхождению и в конечном счете удалению.

Сводные данные

В этом примере вы изучили, как использовать satelliteScenario объект от Aerospace Toolbox™, чтобы импортировать информацию об орбите из файлов TLE. Вы распространили спутниковые траектории с помощью SGP4 и визуализировали сценарий с помощью Спутникового Средства просмотра Сценария. Вы изучили, как использовать модели радара и средства отслеживания от Sensor Fusion and Tracking Toolbox™, чтобы смоделировать радарную систему слежения наблюдения за космическим пространством. Созданная система слежения может предсказать предполагаемую орбиту каждого спутника с помощью модели низкого качества.

Вспомогательные Функции

initKeplerUKF инициализирует Сигма-точечный фильтр Калмана с помощью Кеплеровской модели движения.

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,...
    '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 используется в примере, чтобы определить, какие дорожки обнаруживаемы в установленный срок.

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

deleteBadTracks используется, чтобы удалить дорожки, которые, очевидно, отличались. А именно, отличенные дорожки в этом примере являются дорожками, текущее положение которых упало на поверхность земли и отслеживает, чья ковариация стала слишком большой.

function tracks = deleteBadTracks(tracker,tracks)
% remove divergent tracks:
%   - tracks with covariance > 1e10 (100 km standard deviation)
%   - tracks with estimated position outside of LEO bounds
n = numel(tracks);
toDelete = zeros(1,n,'logical');
for i=1:numel(tracks)
    [pos, cov] = getTrackPositions(tracks(i),[ 1 0 0 0 0 0 ; 0 0 1 0 0 0; 0 0 0 0 1 0]);
    if norm(pos) < 6500*1e3 || norm(pos) > 8500*1e3 || max(cov,[],'all') > 1e10
        deleteTrack(tracker,tracks(i).TrackID);
        toDelete(i) =true;
    end
end
tracks(toDelete) = [];
end

assembleRadarInput используется, чтобы создать положения созвездия в каждой радарной системе координат тела. Это - первый шаг, описанный в схеме.

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 LLA position
%                              and NED 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 convert the current satellite pose to ECEF
targetPoses = repmat(targetTemplate, 1, numel(constellation));
for i=1:numel(constellation)
    lla = constellation(i).Position;
    targetPoses(i).Position = lla2ecef(lla);
    vned = constellation(i).Velocity;
    Rned2ecef = dcmecef2ned(lla(1),lla(2))';
    targetPoses(i).Velocity(:) = Rned2ecef*vned(:);
    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 шаг 2 реализаций в радарном процессе цепи как описано в схеме.

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 используется с метриками присвоения, чтобы оценить присвоение отслеживания.

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

Ссылка

[1] Sridharan, Рамасвами, и Антонио Ф. Пенса, перспективы редакторов в Наблюдении за космическим пространством. Нажатие MIT, 2017.

Для просмотра документации необходимо авторизоваться на сайте