Отследите обломки пробела Используя кеплеровскую модель движения
В этом примере показано, как смоделировать центральные землей траектории с помощью пользовательских моделей движения в trackingScenario, как сконфигурировать fusionRadarSensor в моностатическом режиме, чтобы сгенерировать синтетические обнаружения обломков пробела, и как установить мультиобъектное средство отслеживания, чтобы отследить симулированные цели.
Сценарий обломков пробела
Существует больше чем 30 000 больших объектов обломков (с диаметром, больше, чем 10 см) и больше чем 1 миллион меньших объектов обломков в Низкой околоземной орбите (LEO) [1]. Эти обломки могут быть опасными для деятельности человека на пробеле, повредить операционные спутники и обеспечить время чувствительные и дорогостоящие маневры предотвращения. Когда действие пробела увеличивается, уменьшание и контролирование обломков пробела становятся крайне важными.
Можно использовать Sensor Fusion and Tracking Toolbox™, чтобы смоделировать траектории обломков, сгенерировать синтетические радарные обнаружения этих обломков и получить оценки положения и скорости каждого объекта.
Во-первых, создайте сценарий отслеживания и установите случайный seed для повторяемых результатов.
s = rng; rng(2020); scene = trackingScenario('IsEarthCentered',true, 'InitialAdvance', 'UpdateInterval',... 'StopTime', 3600, 'UpdateRate', 0.1);
Вы используете систему координат Земли фиксированной земли в центре (ECEF). Источник этой системы координат находится в центре Земли и точек оси Z к Северному полюсу. Ось X указывает на пересечение экватора и Гринвичского меридиана. Ось Y завершает предназначенную для правой руки систему. Положения платформы и скорости заданы с помощью Декартовых координат в этой системе координат.
Задайте модель движения обломков
helperMotionTrajectory класс, используемый в этом примере, задает траектории объекта обломков с помощью пользовательской функции модели движения.
Траектории объектов пробела, вращающихся вокруг Земли, могут быть аппроксимированы моделью Keplerian, которая принимает, что Земля является массовым точкой телом, и объекты, движущиеся по кругу вокруг земли, имеют незначительные массы. Эффекты высшего порядка в Наземном поле тяготения и экологических воздействиях не составляются. Поскольку уравнение движения описывается в системе координат ECEF, которая является неинерционной системой координат, Кориолис и центростремительные силы составляются.
Ускоряющий вектор объекта обломков ECEF
,
где стандартный гравитационный параметр Земли, вектор положения объекта обломков ECEF, норма радиус-вектора, и Наземный вектор вращения.
Функциональный 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 'ReferenceRange',2e6,... m 'RangeLimits', [0 2e6],... '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;
Визуализируйте основную истину с trackingGlobeViewer
Вы используете trackingGlobeViewer визуализировать все элементы, заданные в сценарии отслеживания: отдельные обломки возражают и их траектории, радарные вентиляторы, радарные обнаружения и дорожки.
f = uifigure; viewer = trackingGlobeViewer(f, 'Basemap','satellite','ShowDroppedTracks',false); % Add info box on top of the globe viewer infotext = simulationInfoText(0,0,0); infobox = uilabel(f,'Text',infotext,'FontColor',[1 1 1],'FontSize',11,... 'Position',[10 20 300 70],'Visible','on'); % Show radar beams on the globe covcon = coverageConfig(scene); plotCoverage(viewer,covcon, 'ECEF'); while advance(scene) infobox.Text = simulationInfoText(scene.SimulationTime, 0, numDebris); plotPlatform(viewer, debris, 'ECEF'); end %Take a snapshot of the scenario snapshot(viewer);
На виртуальном земном шаре вы видите обломки пробела, представленные белыми точками отдельными запаздывающими траекториями, показанными белыми линиями. Большинство сгенерированных объектов обломков находится на орбитах с высокими наклонными углами близко к 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(viewer); % Reduce history depth for debris viewer.PlatformHistoryDepth = 2; % Show 10-sigma covariance ellipses for better visibility viewer.NumCovarianceSigma = 10; plotCoverage(viewer,covcon, 'ECEF'); % Initialize tracks confTracks = objectTrack.empty(0,1); while advance(scene) plotPlatform(viewer, debris); time = scene.SimulationTime; % Generate detections detections = detect(scene); plotDetection(viewer, detections, 'ECEF'); % Generate and update tracks detectableInput = isDetectable(tracker,time, covcon); if isLocked(tracker) || ~isempty(detections) [confTracks, ~, allTracks,info] = tracker(detections,time,detectableInput); confTracks = deleteBadTracks(tracker,confTracks); plotTrack(viewer, confTracks, 'ECEF'); end infobox.Text = simulationInfoText(time, numel(confTracks), numDebris); % Move camera during simulation and take snapshots switch time case 100 campos(viewer,[90 150 5e6]); camorient(viewer,[0 -65 345]); drawnow im1 = snapshot(viewer); case 270 campos(viewer,[60 -120 2.6e6]); camorient(viewer, [20 -45 20]); drawnow case 380 campos(viewer,[60 -120 2.6e6]); camorient(viewer, [20 -45 20]); drawnow im2 = snapshot(viewer); case 400 % reset campos(viewer,[17.3 -67.2 2.400e7]); camorient(viewer, [360 -90 0]); drawnow case 1500 campos(viewer,[54 2.3 6.09e6]); camorient(viewer, [0 -73 348]); drawnow case 1560 im3 = snapshot(viewer); drawnow end end % Restore random seed rng(s);
imshow(im1);
На первом снимке состояния вы видите объект быть прослеженными как дорожка T1 желтого цвета. Этот объект был только обнаружен дважды, которого было недостаточно, чтобы уменьшать неопределенность в дорожке. Поэтому размер его эллипса ковариации является относительно большим. Можно также наблюдать другую дорожку T2 синего цвета, который несколько раз обнаруживается датчиком. В результате его соответствующий эллипс ковариации намного меньше, поскольку больше обнаружений использовалось, чтобы откорректировать оценку состояния.
imshow(im2);
Несколько минут спустя, как замечено на снимке состояния выше, T1 был удален, потому что неопределенность в дорожке стала слишком большой без обнаружений. С другой стороны, вторая дорожка T2, переживший из-за дополнительных обнаружений.
imshow(im3)
В снимке экрана выше, вы видите, что дорожка T15 (в голубом) собирается ввести радарную область наблюдения. Отследите T11 (в оранжевом), был только обновлен с обнаружением, которое уменьшало неопределенность в его предполагаемом положении и скорости. С радарной настройкой станции, после 30 минут наблюдения, 18 дорожек были инициализированы и подтверждены из 100 объектов обломков. Если вы увеличите время симуляции, радары покроют 360 градусов в области пробела, и в конечном счете больше обломков может быть прослежено. Различные радарные местоположения станции и настройки могли быть исследованы, чтобы увеличить число отслеживаемых объектов.
Сводные данные
В этом примере вы изучили, как задать вашу собственную модель движения, чтобы переместить платформы в сценарий отслеживания и как использовать их, чтобы установить средство отслеживания. Это позволяет вам применить cочетание датчиков и методы отслеживания, предлагаемые в этом тулбоксе более широкой области значений приложений, такие как проблема моделирования и отслеживания обломков пробела в Земле Земля В центре Фиксированная координатная система координат как показано в этом примере.
Вспомогательные Функции
Модель движения, используемая в этом примере, представлена ниже. Состояние является положениями 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
simulationInfoText используется, чтобы отобразить время симуляции, количество обломков и текущее количество дорожек в сценарии.
function text = simulationInfoText(time,numTracks, numDebris) text = vertcat(string(['Elapsed simulation time: ' num2str(round(time/60)) ' (min)']),... string(['Current number of tracks: ' num2str(numTracks)]),... string(['Total number of debris: ' num2str(numDebris)])); drawnow limitrate end
Ссылки
[1] https://www.esa.int/Safety_Security/Space_Debris/Space_debris_by_the_numbers