Управление воздушным движением

В этом примере показано, как сгенерировать сценарий управления воздушным движением, симулируйте радарные обнаружения от радара наблюдения аэропорта (ASR) и сконфигурируйте средство отслеживания глобального самого близкого соседа (GNN), чтобы отследить симулированные цели с помощью радарных обнаружений. Это позволяет вам оценить различные целевые сценарии, радарные требования и настройки средства отслеживания, не нуждаясь в доступе к дорогостоящему самолету или оборудованию. Этот пример покрывает целый синтетический рабочий процесс данных.

Сценарий управления воздушным движением

Симулируйте башню управления воздушным движением (ATC) и движение целей в сценарии как платформы. Симуляция движения платформ в сценарии управляема trackingScenario.

Создайте trackingScenario и добавьте башню ATC в сценарий.

% Create tracking scenario
scenario = trackingScenario;

% Add a stationary platform to model the ATC tower
tower = platform(scenario);

Радар наблюдения аэропорта

Добавьте радар наблюдения аэропорта (ASR) в башню ATC. Типичная башня ATC имеет смонтированные 15 метров радара над землей. Этот радар сканирует механически в азимуте по фиксированной процентной ставке, чтобы предоставить 360 страховых защит степени около башни ATC. Перечислены общие спецификации для ASR:

  • Чувствительность: 0 dBsm 111 км

  • Механическое Сканирование: Азимут только

  • Механическая частота развертки: 12,5 об/мин

  • Электронное сканирование: 'none'

  • Поле зрения: азимут на 1,4 градуса, вертикальное изменение на 10 градусов

  • Разрешение азимута: 1,4 градуса

  • Разрешение области значений: 135 м

Смоделируйте ASR с указанными выше техническими условиями с помощью monostaticRadarSensor.

rpm = 12.5;
fov = [1.4;10];
scanrate = rpm*360/60;  % deg/s
updaterate = scanrate/fov(1); % Hz

radar = monostaticRadarSensor(1, 'Rotator', ...
    'UpdateRate', updaterate, ...           % Hz
    'FieldOfView', fov, ...                 % [az;el] deg
    'MaxMechanicalScanRate', scanrate, ...  % deg/sec
    'AzimuthResolution', fov(1), ...        % deg
    'ReferenceRange', 111e3, ...            % m
    'ReferenceRCS', 0, ...                  % dBsm
    'RangeResolution', 135, ...             % m
    'HasINS', true, ...
    'DetectionCoordinates', 'Scenario');

% Mount radar at the top of the tower
radar.MountingLocation = [0 0 -15];
tower.Sensors = radar;

Наклоните радар так, чтобы он рассмотрел область, начинающуюся в 2 градусах выше горизонта. Для этого включите вертикальное изменение и установите механические пределы сканирования, чтобы охватить поле зрения вертикального изменения радара, начинающееся в 2 градусах выше горизонта. Поскольку trackingScenario использует Северо-восток вниз (NED) координатная система координат, отрицательные вертикальные изменения соответствуют точкам выше горизонта.

% Enable elevation scanning
radar.HasElevation = true;

% Set mechanical elevation scan to begin at 2 degrees above the horizon
elFov = fov(2);
tilt = 2; % deg
radar.MechanicalScanLimits(2,:) = [-fov(2) 0]-tilt; % deg

Установите поле зрения вертикального изменения быть немного больше, чем вертикальное изменение, заполненное пределами сканирования. Это предотвращает растровое сканирование в вертикальном изменении и наклоняет радар, чтобы указать посреди пределов сканирования вертикального изменения.

radar.FieldOfView(2) = elFov+1e-3;

monostaticRadarSensor область значений моделей и вертикальное изменение смещают из-за атмосферного преломления. Эти смещения становятся более явными на более низких высотах и для целей в больших расстояниях. Поскольку индекс изменений преломления (уменьшения) с высотой, радарные сигналы распространяют вдоль изогнутого контура. Это приводит к радарным целям наблюдения на высотах, которые выше, чем их истинная высота и в областях значений вне их области значений угла обзора.

Добавьте три авиалайнера в секторе управления ATC. Один авиалайнер приближается к ATC с большого расстояния, другой отбывает, и третье летит тангенциальный в башню. Смоделируйте движение этих авиалайнеров на 60-секундном интервале.

trackingScenario использует Северо-восток вниз (NED) координатная система координат. При определении waypoints для авиалайнеров ниже, z-координата соответствует вниз, таким образом, высоты над землей установлены в отрицательные величины.

% Duration of scenario
sceneDuration = 60; % s

% Inbound airliner
ht = 3e3;
spd = 900*1e3/3600; % m/s
wp = [-5e3 -40e3 -ht;-5e3 -40e3+spd*sceneDuration -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Outbound airliner
ht = 4e3;
spd = 700*1e3/3600; % m/s
wp = [20e3 10e3 -ht;20e3+spd*sceneDuration 10e3 -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Tangential airliner
ht = 4e3;
spd = 300*1e3/3600; % m/s
wp = [-20e3 -spd*sceneDuration/2 -ht;-20e3 spd*sceneDuration/2 -ht];
traj = waypointTrajectory('Waypoints',wp,'TimeOfArrival',[0 sceneDuration]);
platform(scenario,'Trajectory', traj);

% Create a display to show the true, measured, and tracked positions of the
% airliners.
[theater,fig] = helperATCExample('Create Display',scenario,tower,radar);
helperATCExample('Update Display',theater,scenario,tower,radar);

Средство отслеживания GNN

Создайте trackerGNN сформировать дорожки из радарных обнаружений, сгенерированных от этих трех авиалайнеров. Обновите средство отслеживания с обнаружениями, сгенерированными после завершения полных 360 сканирований степени в азимуте.

Средство отслеживания использует initFilter поддерживание функции, чтобы инициализировать постоянную скорость расширило Фильтр Калмана для каждого нового трека. initFilter изменяет фильтр, возвращенный initcvekf чтобы совпадать с целевыми скоростями и средством отслеживания обновляют интервал.

tracker = trackerGNN( ...
    'Assignment', 'Auction', ...
    'AssignmentThreshold',50, ...
    'FilterInitializationFcn',@initFilter);

Симулируйте и отследите авиалайнеры

Следующий цикл совершенствует положения платформы, пока конец сценария не был достигнут. Для каждого шага вперед в сценарии радар генерирует обнаружения от целей в его поле зрения. Средство отслеживания обновляется с этими обнаружениями после того, как радар завершил 360 сканирований степени в азимуте.

% Set simulation to advance at the update rate of the radar
scenario.UpdateRate = radar.UpdateRate;

% Create a buffer to collect the detections from a full scan of the radar
scanBuffer = {};

% Initialize the track array
tracks = [];

% Set random seed for repeatable results
rng(2018)

while advance(scenario) && ishghandle(fig)

    % Current simulation time
    simTime = scenario.SimulationTime;

    % Target poses in the ATC's coordinate frame
    targets = targetPoses(tower);

    % Use the tower's true position as its INS measurement
    ins = pose(tower, 'true');

    % Generate detections on target's in the radar's current field of view
    [dets,~,config] = radar(targets,ins,simTime);

    scanBuffer = [scanBuffer;dets]; %#ok<AGROW>

    % Update tracks when a 360 degree scan is complete
    if config.IsScanDone
        % Update tracker
        tracks = tracker(scanBuffer,simTime);

        % Clear scan buffer for next scan
        scanBuffer = {};
    elseif isLocked(tracker)
        % Predict tracks to the current simulation time
        tracks = predictTracksToTime(tracker,'confirmed',simTime);
    end

    % Update display with current beam position, buffered detections, and
    % track positions
    helperATCExample('Take Snapshots',fig,scenario,config);
    helperATCExample('Update Display',theater,scenario,tower,radar,scanBuffer,tracks);
    helperATCExample('Take Snapshots',fig,scenario,config);

end
helperATCExample('Take Snapshots',fig,scenario,config);

Используйте helperATCExample показать первый снимок состояния, взятый при завершении второго сканирования радара.

helperATCExample('Show Snapshots',1);

Предыдущий рисунок показывает сценарий в конце вторых 360 сканирований степени радара. Радарные обнаружения, показавшие синими точками, присутствуют для каждого из симулированных авиалайнеров. На данном этапе средство отслеживания было уже обновлено одним полным сканированием радара. Внутренне, средство отслеживания инициализировало дорожки для каждого из авиалайнеров. Эти дорожки покажут после обновления после этого сканирования, когда дорожки будут продвинуты на подтвержденный, удовлетворяя требование подтверждения средства отслеживания 2 хитов из 3 обновлений.

Следующие два снимка состояния показывают отслеживание исходящего авиалайнера.

helperATCExample('Show Snapshots',[2 3]);

Предыдущие рисунки показывают изображение дорожки прежде (слева) и сразу после (справа) обновления средства отслеживания после второго сканирования радара. Обнаружение в фигуре слева используется, чтобы обновить и подтвердить инициализированную дорожку от обнаружения предыдущего сканирования для этого авиалайнера. Рисунок справа показывает положение и скорость подтвержденной дорожки. Предполагаемая скорость дорожки обозначается расширением черной линии от местоположения дорожки. Неопределенность в оценке положения дорожки показывается серым эллипсом. Только после двух обнаружений средство отслеживания установило точную оценку положения и скорости исходящего авиалайнера. Истинная высота авиалайнера составляет 4 км, и она перемещается на восток на уровне 700 км/час.

helperATCExample('Show Snapshots',[4 5]);

Состояние дорожки исходящего авиалайнера курсируется в конец третьего сканирования и показанное на рисунке выше слева наряду с новым обнаружением для авиалайнера. Заметьте, как неопределенность дорожки выросла, поскольку она была обновлена в предыдущей фигуре. Дорожка после того, как это было обновлено с обнаружением, показана на рисунке справа. Вы замечаете, что неопределенность в положении дорожки уменьшается после обновления. Неопределенность дорожки растет между обновлениями и уменьшается каждый раз, когда она обновляется с новым измерением. Вы также замечаете, что после третьего обновления, дорожка находится сверх истинного положения авиалайнера.

helperATCExample('Show Snapshots',6);

Итоговый рисунок показывает состояние дорожек всех трех авиалайнеров в конце сценария. Существует точно одна дорожка для каждого из этих трех авиалайнеров. Те же номера дорожек присвоены каждому из авиалайнеров на целое время сценария, указав, что ни одна из этих дорожек не была пропущена во время сценария. Предполагаемые дорожки тесно совпадают с истинным положением и скоростью авиалайнеров.

truthTrackTable = tabulateData(scenario, tracks) %#ok<NOPTS>
truthTrackTable =

  3x4 table

    TrackID        Altitude              Heading               Speed      
               True    Estimated    True    Estimated    True    Estimated
    _______    _________________    _________________    _________________

     "T1"      4000      4119        90         90       700        701   
     "T2"      4000      4029         0        359       300        285   
     "T3"      3000      3082         0        359       900        890   

Визуализируйте дорожки в 3D, чтобы получить лучший смысл предполагаемых высот.

axis square
view(50,10)
xlim([-25 35])
ylim([-40 20])
zlim([-6 0])

Визуализируйте дорожки на карте

Постройте долготу широты каждого положения дорожки в географические оси с помощью geoplot. Можно использовать ned2geodetic функция от Mapping Toolbox™, чтобы преобразовать положения дорожки от Декартовых координат NED до геодезических координат. Широта и долгота всех дорожек хранятся в ATCExample.mat файл.

helperATCExample('geoPlotTracks','ATCExample.mat');

Сводные данные

В этом примере показано, как сгенерировать сценарий управления воздушным движением, симулируйте радарные обнаружения от радара наблюдения аэропорта (ASR) и сконфигурируйте средство отслеживания глобального самого близкого соседа (GNN), чтобы отследить симулированные цели с помощью радарных обнаружений. В этом примере вы изучили, как основанная на истории логика средства отслеживания продвигает дорожки. Вы также изучили, как неопределенность дорожки растет, когда дорожка курсируется и уменьшается, когда дорожка обновляется новым обнаружением.

Вспомогательные Функции

initFilter

Эта функция изменяет функциональный initcvekf обрабатывать более высокие скоростные цели, такие как авиалайнеры в сценарии ATC.

function filter = initFilter(detection)
filter = initcvekf(detection);
classToUse = class(filter.StateCovariance);

% Airliners can move at speeds around 900 km/h. The velocity is
% initialized to 0, but will need to be able to quickly adapt to
% aircraft moving at these speeds. Use 900 km/h as 1 standard deviation
% for the initialized track's velocity noise.
spd = 900*1e3/3600; % m/s
velCov = cast(spd^2,classToUse);
cov = filter.StateCovariance;
cov(2,2) = velCov;
cov(4,4) = velCov;
filter.StateCovariance = cov;

% Set filter's process noise to match filter's update rate
scaleAccelHorz = cast(1,classToUse);
scaleAccelVert = cast(1,classToUse);
Q = blkdiag(scaleAccelHorz^2, scaleAccelHorz^2, scaleAccelVert^2);
filter.ProcessNoise = Q;
end

tabulateData

Эта функция возвращает таблицу, сравнивающую основную истину и дорожки

function truthTrackTable = tabulateData(scenario, tracks)
% Process truth data
platforms = scenario.Platforms(2:end); % Platform 1 is the radar
numPlats = numel(platforms);
trueAlt = zeros(numPlats,1);
trueSpd = zeros(numPlats,1);
trueHea = zeros(numPlats,1);
for i = 1:numPlats
    traj = platforms{i}.Trajectory;
    waypoints = waypointInfo(traj);
    times = waypoints.TimeOfArrival;
    waypoints = waypoints.Waypoints;
    trueAlt(i) = -waypoints(end,3);
    trueVel = (waypoints(end,:) - waypoints(end-1,:)) / (times(end)-times(end-1));
    trueSpd(i) = norm(trueVel) * 3600 / 1000; % Convert to km/h
    trueHea(i) = atan2d(trueVel(1),trueVel(2));
end
trueHea = mod(trueHea,360);

% Associate tracks with targets
atts = [tracks.ObjectAttributes];
tgtInds = [atts.TargetIndex];

% Process tracks assuming a constant velocity model
numTrks = numel(tracks);
estAlt = zeros(numTrks,1);
estSpd = zeros(numTrks,1);
estHea = zeros(numTrks,1);
truthTrack = zeros(numTrks,7);
for i = 1:numTrks
    estAlt(i) = -round(tracks(i).State(5));
    estSpd(i) = round(norm(tracks(i).State(2:2:6)) * 3600 / 1000); % Convert to km/h;
    estHea(i) = round(atan2d(tracks(i).State(2),tracks(i).State(4)));
    estHea(i) = mod(estHea(i),360);
    platID = tgtInds(i);
    platInd = platID - 1;
    truthTrack(i,:) = [tracks(i).TrackID, trueAlt(platInd), estAlt(i), trueHea(platInd), estHea(i), ...
        trueSpd(platInd), estSpd(i)];
end

% Organize the data in a table
names = {'TrackID','TrueAlt','EstimatedAlt','TrueHea','EstimatedHea','TrueSpd','EstimatedSpd'};
truthTrackTable = array2table(truthTrack,'VariableNames',names);
truthTrackTable = mergevars(truthTrackTable, (6:7), 'NewVariableName', 'Speed', 'MergeAsTable', true);
truthTrackTable.(6).Properties.VariableNames = {'True','Estimated'};
truthTrackTable = mergevars(truthTrackTable, (4:5), 'NewVariableName', 'Heading', 'MergeAsTable', true);
truthTrackTable.(4).Properties.VariableNames = {'True','Estimated'};
truthTrackTable = mergevars(truthTrackTable, (2:3), 'NewVariableName', 'Altitude', 'MergeAsTable', true);
truthTrackTable.(2).Properties.VariableNames = {'True','Estimated'};
truthTrackTable.TrackID = "T" + string(truthTrackTable.TrackID);
end
Для просмотра документации необходимо авторизоваться на сайте