Отслеживайте космический мусор с помощью кеплеровой модели движения

Открыть Live Script

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

Сценарий космического мусора

На околоземной орбите (НОО) насчитывается более 30 000 крупных объектов мусора (диаметром более 10 см) и более 1 млн объектов меньшего размера [1]. Этот мусор может быть опасен для человеческой деятельности в пространстве, повредить эксплуатационные спутники, а также заставить чувствительные и дорогостоящие маневры по избеганию. По мере расширения космической деятельности решающее значение приобретает сокращение и мониторинг космического мусора.

Можно использовать Sensor Fusion and Tracking Toolbox™, чтобы смоделировать траектории обломков, сгенерировать синтетические радиолокационные обнаружения этого мусора и получить оценки положения и скорости каждого объекта.

Сначала создайте сценарий отслеживания и установите случайный seed для повторяемых результатов.

s = rng;
rng(2020);
scene = trackingScenario('IsEarthCentered',true,'InitialAdvance','UpdateInterval');

Вы используете опорную систему координат Earth-Centered-Earth-Fixed (ECEF). Источник этой системы координат находится в центре Земли, а ось Z указывает на северный полюс. Ось X указывает в сторону пересечения экватора и меридиана Гринвича. Ось Y завершает правую систему. Положения и скорости платформы задаются с помощью Декартовых координат в этой системе координат.

Задайте модель движения мусора

The helperMotionTrajectory класс, используемый в этом примере, задает траектории объектов мусора с помощью пользовательской функции модели движения.

Траектории космических объектов, вращающихся вокруг Земли, могут быть аппроксимированы кеплеровской моделью, которая принимает, что Земля является телом точечной массы, а объекты, вращающиеся вокруг Земли, имеют незначительные массы. Эффекты более высокого порядка в гравитационном поле Земли и нарушений порядка окружающей среды не учитываются. Поскольку уравнение движения выражено в системе координат ECEF, которая является неинерционной системой отсчета, учитываются силы Кориолиса и центростремления.

Вектор ускорения объекта мусора ECEF

$\vec{a} = \frac{- μ}{r^{3}} \vec{r} - 2 \vec{Ω} \times \frac{d}{d t} \vec{r} - \vec{Ω} \times (\vec{Ω} \times \vec{r})$ ,

где $μ$ - стандартный гравитационный параметр Земли, $\vec{r}$ - вектор положения объекта мусора ECEF, $r$ является нормой вектора положения, и $\vec{Ω}$ - вектор вращения Земли.

Функция keplerorbit приведенное ниже использование численного интегрирования Рунге-Кутты 4-го порядка этого уравнения для распространения положения и скорости во времени.

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

% Generate a population of debris
numDebris = 100;

range = 7e6 + 1e5*randn(numDebris,1);
ecc = 0.015 + 0.005*randn(numDebris,1);
inc = 80 + 10*rand(numDebris,1);
lan = 360*rand(numDebris,1);
w = 360*rand(numDebris,1);
nu = 360*rand(numDebris,1);

% Convert to initial position and velocity
for i = 1:numDebris
    [r,v] = oe2rv(range(i),ecc(i),inc(i),lan(i),w(i),nu(i));
    data(i).InitialPosition = r; %#ok<SAGROW>
    data(i).InitialVelocity = v; %#ok<SAGROW>
end

% Create platforms and assign them trajectories using the keplerorbit motion model
for i=1:numDebris
    debris(i) = platform(scene); %#ok<SAGROW> 
    debris(i).Trajectory = helperMotionTrajectory(@keplerorbit,...
        'SampleRate',0.1,... % integration step 10sec
        'Position',data(i).InitialPosition,...
        'Velocity',data(i).InitialVelocity); %#ok<SAGROW> 
end

Моделируйте радары наблюдения пространства

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

% Create a space surveillance station in the Pacific ocean
station1 = platform(scene,'Position',[10 180 0]);

% Create a second surveillance station in the Atlantic ocean
station2 = platform(scene,'Position',[0 -20 0]);

% Near the North Pole, create a third surveillance station in Iceland
station3 = platform(scene,'Position',[65 -20 0]);

% Create a fourth surveillance station near the south pole
station4 = platform(scene,'Position',[-90 0 0]);

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

Обнаружение объекта 10 дБсм на расстоянии до 2000 км
Разрешение объектов горизонтально и вертикально с точностью 100 м на области значений 2000 км
Имеет веерообразное поле зрения 120 степеней по азимуту и 30 степени по повышению
Взгляд в пространство на основе его географического положения

% Create fan-shaped monostatic radars to monitor space debris objects
radar1 = fusionRadarSensor(1,...
    'UpdateRate',0.1,... 10 sec
    'ScanMode','No scanning',...
    'MountingAngles',[0 90 0],... look up
    'FieldOfView',[120;30],... degrees
    'RangeLimits', [0 2e6], ... m
    'ReferenceRange',2e6,... m
    'ReferenceRCS', 10,... dBsm
    'HasFalseAlarms',false,...
    'HasElevation',true,...
    'AzimuthResolution',0.01,... degrees
    'ElevationResolution',0.01,... degrees
    'RangeResolution',100,... m
    'HasINS',true,...
    'DetectionCoordinates','Scenario');
station1.Sensors = radar1;

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

radar3 = clone(radar1);
radar3.SensorIndex = 3;
station3.Sensors = radar3;

radar4 = clone(radar1);
radar4.SensorIndex = 4;
station4.Sensors = radar4;

Визуализация основной истины на виртуальном глобусе

В этом примере helperScenarioGlobeViewer предоставляет виртуальный глобус, используемый для визуализации всех элементов, определенных в сценарии слежения: отдельных объектов мусора и их траекторий, радиолокационных вентиляторов, радиолокационных обнаружений и дорожек.

globeDisplay = helperScenarioGlobeViewer;

% Show radar beams on the globe
covcon = coverageConfig(scene);
plotCoverage(globeDisplay,covcon);

% Set TargetHistoryLength to visualize the full trajectory of the debris objects
globeDisplay.TargetHistoryLength = 50;

scene.StopTime = 3600;
scene.UpdateRate = 0.1;

while advance(scene)
    time = scene.SimulationTime;
    updateDisplay(globeDisplay,time,debris);
end
snap(globeDisplay);

На виртуальном глобусе можно увидеть космический мусор, представленный белыми точками с отдельными конечными траекториями, показанными белыми линиями. Большинство сгенерированных объектов мусора находятся на орбитах с высокими углами наклона, близкими к 80 o.

Траектории нанесены в координатах ECEF, и поэтому вся траектория вращается к западу из-за вращения Земли. После нескольких периодов орбиты весь космический мусор проходит через лучи наблюдения радаров.

Моделирование синтетических обнаружений и отслеживание космического мусора

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

Многообъектный трекер trackerJPDA используется для создания новых треков, связи обнаружений с существующими треками, оценки их состояния и удаления расходящихся треков. Установка свойства HasDetectableTrackIDsInput в true позволяет трекеру принять вход, который указывает, является ли отслеживаемый объект обнаруживаемым в области наблюдения. Это важно для того, чтобы не штрафовать дорожки, которые распространяются вне зон радиолокационного наблюдения. Область служебной функции isDetectable вычисляет, какие дорожки обнаруживаются на каждом шаге симуляции.

Кроме того, служебная функция deleteBadTracks используется для быстрого удаления расходящихся дорожек.

% Define Tracker
tracker = trackerJPDA('FilterInitializationFcn',@initKeplerUKF,...
    'HasDetectableTrackIDsInput',true,...
    'ClutterDensity',1e-20,...
    'AssignmentThreshold',1e3,...
    'DeletionThreshold',[7 10]);

% Reset scenario, seed, and globe display
restart(scene);
scene.StopTime = 1800; % 30 min
clear(globeDisplay);
globeDisplay.TargetHistoryLength = 2;
% Scale up size of covariance ellipses for better visibility
globeDisplay.DetectionCovarianceSize = 10;
globeDisplay.TrackCovarianceSize = 10;
plotCoverage(globeDisplay,covcon);

% Initialize tracks
confTracks = objectTrack.empty(0,1);
while advance(scene)
    time = scene.SimulationTime;
    
    % Generate detections
    detections = detect(scene);
    
    % Generate and update tracks
    detectableInput = isDetectable(tracker,time, covcon);
    if ~isempty(detections) || isLocked(tracker)
        [confTracks, ~, allTracks,info] = tracker(detections,time,detectableInput);
        confTracks = deleteBadTracks(tracker,confTracks);
    end

    % Update globe display
    updateDisplay(globeDisplay,time,debris,detections,[],confTracks);
    
    % Move camera during simulation and take snapshots
    switch time
        case 100
            setCamera(globeDisplay,[90 150 5e6],[0 -65 345]);
            im1 = snap(globeDisplay);
        case 270
            setCamera(globeDisplay,[60 -120 2.6e6],[20 -45 20]);
        case 380
            setCamera(globeDisplay,[60 -120 2.6e6],[20 -45 20]);
            im2 = snap(globeDisplay);
        case 400
            % reset
            setCamera(globeDisplay,[17.3 -67.2 2.400e7], [360 -90 0]);
        case 1500
            setCamera(globeDisplay,[54 2.3 6.09e6], [0 -73 348]);
        case 1560
            im3 = snap(globeDisplay);
    end
        
end

% Restore random seed
rng(s);

imshow(im1);

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

imshow(im2);

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

imshow(im3)

На скриншоте выше можно увидеть, что трек T15 (желтым цветом) собирается войти в зону радиолокационного наблюдения. Трек- T11 (в зеленом цвете) был только что обновлен обнаружением, что уменьшило неопределенность его предполагаемого положения и скорости. С строением радиолокационной станции после 30 минут наблюдения 18 пути были инициализированы и подтверждены из 100 объектов обломков. Если вы увеличите время симуляции, радары накроют 360 степени в пространстве и в конечном счете можно будет отслеживать больше обломков. Различные местоположения и строения радиолокационных станций могут быть исследованы, чтобы увеличить количество отслеживаемых объектов.

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

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

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

Модель движения, используемая в этом примере, представлена ниже. Состояние является положениями ECEF и скоростями объекта [x; vx; y; vy; z; vz].

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 4 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

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

function filter = initKeplerUKF(detection)

% assumes radar returns [x y z]
measmodel= @(x,varargin) x([1 3 5],:);
detNoise = detection.MeasurementNoise;
sigpos = 0.4;% m
sigvel = 0.4;% m/s^2
meas = detection.Measurement;
initState = [meas(1); 0; meas(2); 0; meas(3);0];
filter = trackingUKF(@keplerorbit,measmodel,initState,...
    'StateCovariance', diag(repmat([10, 10000].^2,1,3)),...
    'ProcessNoise',diag(repmat([sigpos, sigvel].^2,1,3)),...
    'MeasurementNoise',detNoise);

end

oe2rv преобразует набор из 6 традиционных орбитальных элементов в векторы положения и скорости.

function [r,v] = oe2rv(a,e,i,lan,w,nu)
% Reference: Bate, Mueller & White, Fundamentals of Astrodynamics Sec 2.5

mu = 398600.4405*1e9; % m^3 s^-2

% Express r and v in perifocal system
cnu = cosd(nu);
snu = sind(nu);
p = a*(1 - e^2);
r = p/(1 + e*cnu);
r_peri = [r*cnu ; r*snu ; 0];
v_peri = sqrt(mu/p)*[-snu ; e + cnu ; 0];

% Tranform into Geocentric Equatorial frame
clan = cosd(lan);
slan = sind(lan);
cw = cosd(w);
sw = sind(w);
ci = cosd(i);
si = sind(i);
R = [ clan*cw-slan*sw*ci  ,  -clan*sw-slan*cw*ci   ,    slan*si; ...
    slan*cw+clan*sw*ci  ,   -slan*sw+clan*cw*ci  ,   -clan*si; ...
    sw*si                  ,   cw*si                   ,   ci];
r = R*r_peri;
v = R*v_peri;
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 > 4*1e8 (20 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') > 4*1e8
        deleteTrack(tracker,tracks(i).TrackID);
        toDelete(i) =true;
    end
end
tracks(toDelete) = [];
end

Ссылки

[1] https://www.esa.int/Safety_Security/Space_Debris/Space_debris_by_the_numbers

Документация