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