Моделирование и отслеживание маршрутных самолетов в сценариях, ориентированных на Землю

В этом примере показано, как использовать центрированную землей 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

Ссылка

[1] FAR § 91.227 Автоматическое зависимое наблюдение-широковещание (ADS-B) Требования к эффективности оборудования