В этом примере показано, как смоделировать центральные землей траектории с помощью пользовательских моделей движения в 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 визуализировать все элементы, заданные в сценарии отслеживания: отдельные обломки возражают и их траектории, радарные вентиляторы, радарные обнаружения и дорожки.
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
