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

Скрипт Open Live Script

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

Сценарий обломков пробела

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

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

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

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

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

Задайте модель движения обломков

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

Траектории объектов пробела, вращающихся вокруг Земли, могут быть аппроксимированы моделью Keplerian, которая принимает, что Земля является массовым точкой телом, и объекты, движущиеся по кругу вокруг земли, имеют незначительные массы. Эффекты высшего порядка в Наземном поле тяготения и экологических воздействиях не составляются. Поскольку уравнение движения описывается в системе координат 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 объектов dBsm на расстоянии в 2 000 км
Решение объектов горизонтально и вертикально с точностью 100 м в 2 000-километровой области значений
Наличие веерообразного поля зрения 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 градусам.

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

Симулируйте синтетические обнаружения и отследите обломки пробела

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

Мультиобъектное средство отслеживания trackerJPDA используется, чтобы создать новые треки, объединенные обнаружения к существующим дорожкам, оценить их состояние и удалить расходящиеся дорожки. Установка свойства HasDetectableTrackIDsInput к истине позволяет средству отслеживания принимать вход, который указывает, обнаруживаем ли отслеживаемый объект в области наблюдения. Это важно для того, чтобы не оштрафовать дорожки, которые распространены за пределами радарных областей наблюдения. Служебная функция 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

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