В этом примере показано, как моделировать траектории, ориентированные на землю, с использованием пользовательских моделей движения в пределах trackingScenario, как настроить fusionRadarSensor в моностатическом режиме для создания синтетических обнаружений космического мусора и как установить многообъектный трекер для отслеживания имитируемых целей.
На низкой околоземной орбите (НОО) насчитывается более 30 000 крупных объектов мусора (диаметром более 10 см) и более 1 миллиона более мелких объектов мусора [1]. Этот мусор может быть опасен для деятельности человека в космосе, повредить эксплуатационные спутники, а также заставить время чувствительные и дорогостоящие маневры избегания. По мере роста космической деятельности решающее значение приобретает сокращение космического мусора и контроль за ним.
Вы можете использовать Sensor Fusion and Tracking Toolbox™ для моделирования траекторий обломков, создания синтетических радиолокационных обнаружений этого мусора и получения оценок положения и скорости каждого объекта.
Сначала создайте сценарий отслеживания и задайте случайное начальное значение для воспроизводимых результатов.
s = rng; rng(2020); scene = trackingScenario('IsEarthCentered',true,'InitialAdvance','UpdateInterval');
Используется опорная рамка Earth-Centered-Earth-Fixed (ECEF). Начало этого кадра находится в центре Земли, а ось Z указывает на северный полюс. Ось X указывает на пересечение экватора и Гринвичского меридиана. Ось Y завершает правую систему. Положения и скорости платформы определяются с использованием декартовых координат в этом кадре.
helperMotionTrajectory класс, используемый в этом примере, определяет траектории объектов мусора с помощью пользовательской функции модели движения.
Траектории космических объектов, вращающихся вокруг Земли, можно аппроксимировать кеплеровской моделью, которая предполагает, что Земля является телом точечной массы, а объекты, вращающиеся вокруг Земли, имеют незначительные массы. Эффекты более высокого порядка в гравитационном поле Земли и возмущениях окружающей среды не учитываются. Поскольку уравнение движения выражается в кадре ECEF, который является неинерциальным опорным кадром, учитываются кориолисовые и центростремительные силы.
Вектор ускорения объекта мусора ECEF:
),
где - стандартный гравитационный параметр Земли, - вектор положения объекта обломков ECEF, - норма вектора положения, а Ω→is вектор вращения Земли.
Функция 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 объект. Для обнаружения объектов мусора в диапазоне НОО РЛС предъявляет следующие требования:
Обнаружение объекта 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 °.
Траектории строятся в координатах 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