Моделирование и отслеживание маршрутных самолетов в сценариях, ориентированных на Землю
В этом примере показано, как использовать центрированную землей trackingScenario и a geoTrajectory объект для моделирования траектории рейса, которая охватывает тысячи километров. Вы используете две различные модели, чтобы сгенерировать синтетические обнаружения самолета: моностатический радар и отчеты ADS-B. Вы используете трекер мультиобъекта, чтобы оценить траекторию плоскости, сравнить эффективность отслеживания и исследовать влияние, которое ADS-B оказывает на общее качество отслеживания.
В Соединенных Штатах Федеральное управление авиации (FAA) отвечает за регулирование тысяч рейсов каждый день через все национальное воздушное пространство. Коммерческие рейсы обычно отслеживаются в любое время, от их аэропорта вылета до прибытия. Система управления воздушным движением является комплексной многоуровневой системой. Диспетчерские вышки аэропорта отвечают за мониторинг и обработку в непосредственной близости от аэропорта, в то время как центры управления движением воздушных маршрутов (ARTCC) несут ответственность за длительную область значений оперативное наблюдение за различными областями, которые составляют национальное воздушное пространство.
Возможности, а также сложность радаров воздушного сообщения/наблюдения значительно возросли за последние десятилетия. Сложение транспондеров на самолете добавляет двухстороннюю связь между радиолокационным средством и самолетами, что позволяет очень точно оценить положение и выгоды принятия решений в центрах управления. В 2020 году все самолеты, летающие выше 10 000 футов, должны быть оснащены транспондером автоматического зависимого видеотрансляции (ADS-B) для трансляции их бортового предполагаемого положения. Это сообщение принимается и обрабатываются контроллеры воздушного движения.
Создайте сценарий маршрутного воздушного движения
Вы начинаете с создания сценария, ориентированного на Землю.
% Save current random generator state s = rng; % Set random seed for predictable results rng(2020); % Create scenario scene = trackingScenario('IsEarthCentered',true,'UpdateRate',1);
Сгенерируйте истинную траекторию самолета
Матрица flightwaypoints прилагаемый в этом примере содержит синтезированные координаты и временную информацию траектории рейса из Уичиты в Чикаго. Вы используете geoTrajectory объект для создания траектории рейса.
load('flightwaypoints.mat') flightRoute = geoTrajectory(lla,timeofarrival); airplane = platform(scene,'Trajectory',flightRoute);
В настоящее время все коммерческие самолеты оснащены GPS- приемников в составе своих бортовых приборов. Чтобы смоделировать полетные приборы, вы устанавливаете PoseEstimator свойство платформы самолета. В качестве основы ADS-B точность встроенного GPS может быть установлена в соответствии с требованиями ADS-B. Категория навигационной точности, используемая в ADS-B, для положения и скорости называются NACp и NACv. В соответствии с правилами FAA [1] NACp должен быть менее 0,05 морских миль, а NACv - менее 10 метров в секунду. В этом примере вы используете insSensor модель с точностью положения 100 м и точностью скорости 10 м/с.
onboardINS = insSensor('PositionAccuracy',100,'VelocityAccuracy',10,... "RandomStream","mt19937ar with seed"); airplane.PoseEstimator = onboardINS; load('737rcs.mat'); airplane.Signatures{1} = boeing737rcs;
Добавить станции наблюдения вдоль маршрута
Существует несколько моделей радаров дальнего наблюдения, используемых FAA. Радар наблюдения воздушного маршрута 4 (ARSR-4) - радар, введенный в 1990-х годах, который может обеспечить 3D возвраты любого объекта площадью 1 квадратный метр при длине области значений 250 морских миль (463 километра). Большинство ARSR-4 радаров расположены вдоль границ континентальных Соединенных Штатов, в то время как немного более короткие радары дальности в основном расположены на радиолокационных площадках FAA на континенте. В этом примере один тип радара моделируется следующими общими спецификациями ARSR-4, как описано ниже:
Частота обновления: 12 секунд
Максимальная область значений (1 метр квадратная цель): 463 км
Разрешение в области значений: 323 m
Точность области значений: 116 m
Поле зрения азимута: 360 град
Азимут Разрешение: 1.4 град
Точность азимута: 0.176 град
Точность высоты: 900 м
Платформа добавляется к сценарию для каждого радиолокационного узла. Сигнатура RCS этих платформ устанавливается равной -50 децибел, чтобы избежать создания нежелательных радиолокационных возвратов.
По умолчанию радиолокационные обнаружения сообщаются в каркасе кузова радиолокационной монтажной платформы, который в этом случае является локальной системой координат NED в положении каждой радиолокационной площадки. Однако в этом примере вы устанавливаете DetectionCoordinates свойство к Scenario в порядок вывода обнаружений в систему координат ECEF, что позволяет трекеру обрабатывать все обнаружения от различных радиолокационных сайтов в общей системе координат.
% Model an ARSR-4 radar updaterate = 1/12; fov = [360;30.001]; elAccuracy = atan2d(0.9,463); % 900m accuracy @ max range elBiasFraction = 0.1; arsr4 = monostaticRadarSensor(1,'UpdateRate',updaterate,... 'FieldOfView',fov,..., 'HasElevation',true,... 'ScanMode','Mechanical',... 'MechanicalScanLimits',[0 360; -30 0],... 'HasINS',true,... 'HasRangeRate',true,... 'HasFalseAlarms',false,... 'ReferenceRange',463000,... 'ReferenceRCS',0,... 'AzimuthResolution',1.4,... 'AzimuthBiasFraction',0.176/1.4,... 'ElevationResolution',elAccuracy/elBiasFraction,... 'ElevationBiasFraction',elBiasFraction,... 'RangeResolution', 323,... 'RangeBiasFraction',116/323,... Accuracy / Resolution 'RangeRateResolution',100,... 'DetectionCoordinates','Scenario'); % Add ARSR-4 radars at each ARTCC site radarsitesLLA = [41.4228 -88.0583 0;... 40.6989 -89.8258 0;... 39.2219 -95.2461 0]; for i=1:3 radar = clone(arsr4); radar.SensorIndex = i; platform(scene,'Position',radarsitesLLA(i,:),... 'Signatures',rcsSignature('Pattern',-50),'Sensors',{radar}); end
Определите центральный трекер
Как правило, один ARTCC поддерживает дорожки для всех объектов в пределах своей зоны наблюдения и передает информацию отслеживания в следующий ARTCC, когда объект летит в новую зону. В этом примере вы задаете один централизованный трекер для всех радаров. Вы используете trackerGNN объект для сплавления радиолокационных обнаружений плоскости с нескольких радаров с другой информацией о датчике, такой как ADS-B.
tracker = trackerGNN('FilterInitializationFcn',@initfilter,... 'ConfirmationThreshold',[3 5],... 'DeletionThreshold',[5 5],... 'AssignmentThreshold',[1000 Inf]);
Визуализация сцены
Вы используете helperGlobeViewer объект, приложенный в этом примере для отображения платформ, траекторий, обнаружений и дорожек на Земле.
gl = helperGlobeViewer; setCamera(gl,[28.9176 -95.3388 5.8e5],[0 -30 10]); % Show radar sites plotPlatform(gl,scene.Platforms(2:end),'d'); % Show radar coverage covcon = coverageConfig(scene); plotCoverage(gl,covcon); % Show flight route plotTrajectory(gl,airplane); % Take a snapshot snap(gl);
Отслеживайте рейс только с помощью радаров
В этой первой симуляции вы будете симулировать сцену и отслеживать самолет, основанный исключительно на радарах большой области значений.
useADSB = false; snapTimes = [300 1600]; % Declare loop variables detBuffer = {}; tLLA = []; tHeading = []; tSpeed = []; tTime = []; images = {}; % Set updata rates in seconds trackupdatetime = 12; adsbupdatetime = 1; arsrupdatetime = 12; while advance(scene) time = scene.SimulationTime; % Update position of airplane on the globe plotTarget(gl,'airplane',airplane,'^','full'); % Generate radar detections at the defined rate if mod(time,arsrupdatetime) == 0 % Generate synthetic radar detections dets = detect(scene); dets = removeBelowGround(dets); detBuffer = [detBuffer; dets]; %#ok<AGROW> end % Generate ADS-B detections at the defined rate if useADSB && mod(time,adsbupdatetime) == 0 adsb = adsbDetect(time,airplane); detBuffer = [detBuffer; adsb]; %#ok<AGROW> end % Fuse detections in tracker and update tracks on the globe if mod(time,trackupdatetime) == 0 && time > 0 % Update detections on the globe plotDetection(gl,detBuffer); % Tracker needs detections for first call i.e. cond = ~isempty(detBuffer) || ~isLocked(tracker); if cond tracks = tracker(detBuffer,time); if ~isempty(tracks) % Record the estimated airplane data [tLLA, tSpeed, tHeading, tTime] = ... helperNavigationData(tracks,tLLA, tSpeed, tHeading, tTime); end plotTrack(gl,tracks); detBuffer = {}; end end % Move camera and take snapshots images = moveCamera(gl,time,snapTimes,images); end
Обратите внимание, что обнаружения ARSR относительно неточны по высоте, что в целом приемлемо, так как воздушное сообщение контроллера отдельные самолеты горизонтально, а не вертикально. Наименее точные обнаружения генерируются радиолокационными площадками, расположенными на больших областях значений. Эти обнаружения, тем не менее, связаны с дорожкой. На рисунке белая линия представляет истинную траекторию, а желтая - траекторию, оцененную трекером.
for i=1:numel(images) figure imshow(images{i}); end
Первая ветвь рейса далека от любого радара и обнаружения имеют высокий измерительный шум. Кроме того, модель движения с постоянной скоростью не моделирует движение хорошо во время начального поворота рейса после взлета. Когда самолет входит в зону покрытия второго радара наблюдения, дополнительные обнаружения корректируют дорожку, показанную на втором снимке. Трек (желтым цветом) следует истине (белым цветом) намного близко, когда входит в зону покрытия второго радара.
Отслеживайте рейс с помощью радаров и ADS-B
Область служебной функции adsbDetect использует расчетное положение, измеренное собственным навигационным прибором самолета, смоделированным PoseEstimator свойство платформы. Предполагаемое положение обычно кодируется и транслируется на канале 1090 МГц для ближайших приемников ADS-B, чтобы забрать. в США почти полный охват коммерческих рейсов, курсирующих на высотах 10 000 метров и ниже.
useADSB = true; snapTimes = [200 1244]; [trackLLA2, trackSpeed2, trackHeading2, trackTime2, images] = ... simulateScene(gl,scene,tracker,useADSB,snapTimes); % Reset random seed rng(s);
Отображение снимков во время симуляции показано ниже, как и в предыдущем разделе.
for i=1:numel(images) figure imshow(images{i}); end
Обнаружения ADS-B представлены фиолетовым цветом на глобусе. Трасса более тесно совпадает с основной истиной, хотя значительного улучшения при входе в среднюю зону радиолокационного покрытия не наблюдается.
Результаты
Вы сравниваете записанные данные дорожки в обеих симуляциях. Истинное положение и скорость доступны из geoTrajectory. В управлении воздушным движением динамическими состояниями интересов являются положения, представленные широтой, долготой и высотой, курсом и глубиной грунта.
% Query truth at the recorded timestamps for both runs [trueLLA, ~, trueVel] = lookupPose(flightRoute,tTime); [trueLLAadsb, ~, trueVeladsb] = lookupPose(flightRoute,trackTime2); figure tiledlayout(3,1); nexttile hold on plot(tTime/60,trueLLA(:,3),'LineWidth',2,'LineStyle','--','DisplayName','Truth'); plot(tTime/60,tLLA(:,3),'DisplayName','Surveillance radar only'); plot(trackTime2/60,trackLLA2(:,3),'DisplayName','Surveillance radar + ADS-B'); legend('Location','northeastoutside') xlabel('Time (minute)'); ylabel('Altitude (meter)') title('Altitude') nexttile hold on plot(tTime/60,vecnorm(trueVel(:,(1:2)),2,2),'LineWidth',2,'LineStyle','--'); plot(tTime/60,tSpeed); plot(trackTime2/60,trackSpeed2); xlabel('Time (minute)'); ylabel('Speed (m/s)') title('Groundspeed') nexttile hold on plot(tTime/60,atan2d(trueVel(:,2),trueVel(:,1)),'LineWidth',2,'LineStyle','--'); plot(tTime/60,tHeading); plot(trackTime2/60,trackHeading2); xlabel('Time (minute)'); ylabel('Heading (degree)') title('Heading')
Результаты показывают, что дополнительные данные, предоставленные ADS-B, оказывают значительный эффект на оценку высоты, поскольку высота является наименее точной информацией, представленной радарами. Между тем, оценка грунтовой скорости и курса также улучшается с учетом более частых обновлений (каждые секунды по сравнению с каждые 12 секунд).
Сводные данные
В этом примере вы научились создавать земной сценарий и определять траектории с помощью геодезических координат. Вы также научились моделировать радары наблюдения воздушного маршрута и генерировать синтетические обнаружения. Вы подаете эти обнаружения на трекер мультиобъекта и оцениваете положение, скорость и курс плоскости. Эффективность отслеживания повышается за счет добавления и слияния информации ADS-B. Вы использовали простой подход к моделированию отчетов ADS-B и интеграции их с решением для отслеживания. В этом примере был смоделирован только один рейс. Однако преимущество ADS-B может быть дополнительно проявлено при моделировании нескольких рейсов в сценарии, где безопасные расстояния разделения должны поддерживаться контроллерами воздушного движения.
Вспомогательные функции
simulateScene содержит цикл симуляции, который можно запускать с ADS-B или без. Сценарий обновляется каждую секунду. Обнаружение ADS-B генерируется каждую секунду, и обнаружение ARSR доступно в конце каждые 12 секунд скана. Функция выводит расчетное положение плоскости как широта (град), долгота (град) и WGS84 высота (м), расчетную наземную скорость (м/с), расчетный курс (град) и соответствующие времена (времена) симуляции (ы).
function [tLLA,tSpeed,tHeading,tTime,images] = simulateScene(gl,scene,tracker,useADSB,snapTimes) % Reset the scenario, globe, and seed rng(2020); restart(scene); airplane = scene.Platforms{1}; release(tracker); clear(gl); setCamera(gl,[39.4563 -90.1187 2.356e6]); % Display initial state of the scene covcon = coverageConfig(scene); plotPlatform(gl,scene.Platforms(2:end),'d'); plotCoverage(gl,covcon); % Declare loop variables detBuffer = {}; tLLA = []; tHeading = []; tSpeed = []; tTime = []; images = {}; % Set update rates in seconds if useADSB trackupdatetime = 1; else trackupdatetime = 12; end adsbupdatetime = 1; arsrupdatetime = 12; while advance(scene) time = scene.SimulationTime; % Update position of airplane on the globe plotTarget(gl,'airplane',airplane,'^','full'); % Generate radar detections at the defined rate if mod(time,arsrupdatetime) == 0 && time > 0 % Generate synthetic radar detections dets = detect(scene); dets = removeBelowGround(dets); detBuffer = [detBuffer; dets]; %#ok<AGROW> end % Generate ADS-B detections at the defined rate if useADSB && mod(time,adsbupdatetime) == 0 adsb = adsbDetect(time,airplane); detBuffer = [detBuffer; adsb]; %#ok<AGROW> end % Update detections on the globe plotDetection(gl,detBuffer); % Fuse detections in tracker and update tracks on the globe if mod(time,trackupdatetime) == 0 && time > 0 % Tracker needs detections for first call cond = ~(isempty(detBuffer) && ~isLocked(tracker)); if cond tracks = tracker(detBuffer,time); if ~isempty(tracks) % Record the estimated airplane data [tLLA, tSpeed, tHeading, tTime] = ... helperNavigationData(tracks,tLLA, tSpeed, tHeading, tTime); end plotTrack(gl,tracks); detBuffer = {}; end end % Move camera and take snapshots images = moveCamera(gl,time,snapTimes,images); end end
adsbDetect генерирует обнаружение самолета на основе его бортового набора датчиков, смоделированного insSensor.
function detection = adsbDetect(time,airplane,~) %adsbDetect generate ADS-B detections estimpose = pose(airplane,'CoordinateSystem','Cartesian'); meas = [estimpose.Position(:) ; estimpose.Velocity(:)]; measnoise = diag([100 100 100 10 10 10].^2); mparams = struct('Frame','Rectangular',... 'OriginPositin',zeros(3,1),... 'OriginVelocity',zeros(3,1),... 'Orientation',eye(3),... ADS-B does not report orientation 'IsParentToChild', 1,... 'HasAzimuth', 1,... 'HasElevation', 1,... 'HasRange', 1, ... 'HasVelocity', 1); attr ={ struct('TargetIndex',airplane.PlatformID,'SNR',10)}; detection = {objectDetection(time, meas, 'SensorIndex',20,... 'MeasurementNoise',measnoise,... 'MeasurementParameters',mparams,... 'ObjectAttributes',attr)}; end
initfilter определяет расширенный фильтр калмана, используемый trackerGNN. Движение самолета хорошо аппроксимируется моделью движения с постоянной скоростью. Поэтому довольно небольшой технологический шум придаст больше веса динамике по сравнению с измерениями, которые, как ожидается, будут довольно шумными на больших областях значений.
function filter = initfilter(detection) filter = initcvekf(detection); filter.StateCovariance = 4*filter.StateCovariance; % initcvekf uses measurement noise as the default filter.ProcessNoise = 0.02*eye(3); end
removeBelowGround удаляет радиолокационные обнаружения, которые лежат ниже поверхности Земли, так как они не могут быть созданы реальным радаром.
function detsout = removeBelowGround(detsin) n = numel(detsin); keep = zeros(1,n,'logical'); for i=1:n meas = detsin{i}.Measurement(1:3)'; lla = fusion.internal.frames.ecef2lla(meas); if lla(3)>0 keep(i) = true; else keep(i) = false; end end detsout = detsin(keep); end
moveCamera задает новые положения камеры, чтобы следовать за самолетом и делать снимки.
function images = moveCamera(gl, time, snapTimes,images) if time == 120 setCamera(gl,[37, -97.935, 4.328e4],[0 -23 30]); end if time == 1244 setCamera(gl,[38.8693 -95.6109 1.214e4],[0 -26.8 38.7]); end if time == 1445 setCamera(gl,[37.995, -94.60 6.87e4],[0 -15 2]); end if time == 2100 setCamera(gl, [38.9747 -94.6385 1.3e5], [0 -20 55]); end if time == 3000 setCamera(gl,[41.636 -87.011 6.33e4],[0 -25 270]); end % Snaps if any(time == snapTimes) img = snap(gl); images = [ images, {img}]; end end