Создайте в пути сценарий воздушного движения

Вы запускаете путем создания сосредоточенного Землей сценария.

% Save current random generator state
s = rng;
% Set random seed for predictable results
rng(2020);
% Create scenario
scene = trackingScenario('IsEarthCentered',true,'UpdateRate',1);

Задайте модель самолета и траекторию

matfile flightwaypoints присоединенный в этом примере содержит синтезируемые координаты и информацию времени траектории рейса от Уичито до Чикаго. Вы используете geoTrajectory объект создать траекторию рейса.

load('flightwaypoints.mat')
flightRoute = geoTrajectory(lla,timeofarrival);
airplane = platform(scene,'Trajectory',flightRoute);

В наше время коммерческие плоскости все оборудованы GPS-приемниками. Как магистраль ADS-B, точность встроенного GPS может собираться соответствовать требованиям ADS-B. Категория Точности Навигации, используемая в ADS-B, для положения и скорости, упоминается как NACp и NACv, соответственно. На постановление [1] ФАА NACp должен быть меньше 0,05 морских миль, и NACv должен быть меньше 10 метров в секунду. В этом примере вы используете gpsSensor модель с точностью положения 50 м и скоростной точностью 10 м/с, чтобы сконфигурировать adsbTransponder модель. Вы также используете более реалистическую подпись ЭПР для самолета, вдохновленного тем из Boeing 737.

posAccuracy = 50; % meters
velAccuracy = 10; % m/s

gps = gpsSensor('PositionInputFormat','Geodetic','HorizontalPositionAccuracy',posAccuracy,...
    'VerticalPositionAccuracy', posAccuracy,'VelocityAccuracy',10);

adsbTx = adsbTransponder('MW2020','UpdateRate', 1,'GPS', gps,'Category',adsbCategory.Large);

load('737rcs.mat');
airplane.Signatures{1} = boeing737rcs;

Добавьте станции наблюдения вдоль маршрута

Существует несколько моделей радаров наблюдения дальних, используемых ФАА. Радар Наблюдения Авиалинии 4 (ARSR-4) является радаром, введенным в 1990-х, который может обеспечить 3D возвраты любого объекта на 1 квадратный метр в большом расстоянии 250 морских миль (463 километра). Большинство радаров ARSR-4 расположено вдоль границ континентальной части США, в то время как немного короче располагаются, радары в основном расположены на радарных сайтах ФАА на континенте. В этом примере один радарный тип моделируется после общих технических требований ARSR-4, как описано ниже:

Платформа добавляется к сценарию для каждого радарного сайта. Подпись ЭПР тех платформ собирается быть $-$50 децибелов, чтобы не создавать нежелательный радар возвращаются.

По умолчанию о радарных обнаружениях сообщают в радаре, монтирующем систему координат тела платформы, которая в этом случае является локальным Северо-востоком, Вниз структурируют в положении каждого радарного сайта. Однако в этом примере вы устанавливаете 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 = fusionRadarSensor(1,'UpdateRate',updaterate,...
    'FieldOfView',[360 ; 30.001],...
    'HasElevation',true,...
    'ScanMode','Mechanical',...
    'MechanicalAzimuthLimits',[0 360],...
    'MechanicalElevationLimits',[-30 0],...
    'HasINS',true,...
    'HasRangeRate',true,...
    'HasFalseAlarms',false,...
    'ReferenceRange',463000,...
    'ReferenceRCS',0,...
    'RangeLimits', [0 463000],...
    '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

Вы используете adsbReceiver смоделировать прием сообщений ADS-B. Сообщение ADS-B содержит положение, измеренное самолетом собственный инструмент GPS. Сообщение обычно кодируется и широковещательно передается на канале на 1 090 МГц для соседних приемников ADS-B, чтобы взять. Вы задаете область значений приема 200 км вокруг каждой станции наблюдения. В этом примере вы принимаете, что станции наблюдения имеют совершенную связь друг между другом. Поэтому центральный приемник берет, широковещательно передал сообщения ADS-B в области значений по крайней мере одной станции.

% Define central adsbReceiver
adsbRx = adsbReceiver('ReceiverIndex',2);
adsbRange = 2e5;

Визуализируйте сцену

Вы используете trackingGlobeViewer к подиумам, траекториям, обнаружениям и дорожкам на Земле.

viewer = trackingGlobeViewer('Basemap','streets-dark',...
    'TrackLabelScale',1.3,'TrackHistoryDepth',4000,...
    'CoverageMode','Coverage');
% Show radar sites
plotPlatform(viewer,scene.Platforms(2:end));
% Show radar coverage
covcon = coverageConfig(scene);
plotCoverage(viewer,covcon,'ECEF','Alpha',0.1);
% Show flight route
plotTrajectory(viewer,flightRoute);

Радары наблюдения имеют антенну слепой конус, иногда называемый "конусом тишины". Это - объем пробела, непосредственно выше радара, за которым нельзя следить из-за ограничений сканирования антенны. Перекрывающиеся покрытия в сетях радаров являются общей стратегией смягчения этой слепой конической области. С перекрывающейся стратегией, однако, могут все еще быть области, не полностью покрытые сетью. В этом примере конус самого южного радара тишины (отображенный оранжевым в изображении выше) только частично покрыт смежным радаром в сети. Это создает мертвую точку, где самолет не будет обнаружен никаким радаром.

Задайте центральное радарное средство отслеживания и термофиксатор дорожки.

Как правило, один ARTCC обеспечивает дорожки для всех объектов в его зоне наблюдения и передает информацию об отслеживании следующему ARTCC, когда объект летит в новую зону. В этом примере вы задаете одно централизованное средство отслеживания для всех радаров. Вы используете trackerGNN возразите, чтобы плавить радарные обнаружения плоскости от нескольких радаров.

radarGNN = trackerGNN('TrackerIndex',1,...
    'MaxNumTracks',50,...
    'FilterInitializationFcn',@initfilter,...
    'ConfirmationThreshold',[3 5],...
    'DeletionThreshold',[5 5],...
    'AssignmentThreshold',[250 2000]);

Вы также плавите радарные дорожки с дорожками ADS-B, полученными из приемника ADS-B. Для этого вы конфигурируете центральный trackFuser объект. Вы устанавливаете ConfirmationThreshold и DeletionThreshold с учетом различия в частоте обновления между уровнем приемника ADS-B на уровне 1 Гц и радарным уровнем средства отслеживания на уровне 1/12 Гц. Чтобы позволить, чтобы по крайней мере две радарных дорожки были присвоены центральной дорожке, количество хитов должно затем быть, по крайней мере, $2\times 12$.

fuser = trackFuser('FuserIndex',3,'MaxNumSources',2,...
    "AssignmentThreshold",[250 500],...
    "StateFusion",'Intersection',...
    'StateFusionParameters','trace',...
    'ConfirmationThreshold',[2 3*12],...
    'DeletionThreshold',[4 4]*12);

Отследите рейс с помощью радаров и ADS-B

В этом разделе вы симулируете сценарий и шаг радарное средство отслеживания, отслеживаете термофиксатор, транспондер ADS-B и приемник.

% Declare some variables
detBuffer = {};
radarTrackLog = {};
fusedTrackLog = {};
adsbTrackLog = {};
images = {};
snapTimes = [84,564,1128, 2083];

% Track labels and colors
adsblabel = "       ADS-B";
radarlabel = "  Radar";
fusedlabel = string(sprintf('%s\n',"","Fused"));
adsbclr = [183 70 255]/255;
radarclr = [255 255 17]/255;
fusedclr = [255 105 41]/255;

while advance(scene)
    time = scene.SimulationTime;

    % Update position of airplane on the globe
    plotPlatform(viewer,airplane,'TrajectoryMode','None');

    % Generate and plot synthetic radar detections
    [dets, configs] = detect(scene);
    dets = postProcessDetection(dets);
    detBuffer = [detBuffer; dets]; %#ok<AGROW>
    plotDetection(viewer,detBuffer,'ECEF');

    % tracks
    adsbTracks = objectTrack.empty;
    radarTracks = objectTrack.empty;
    fusedTracks = objectTrack.empty;

    % Generate ADS-B messages
    truePose = pose(airplane, 'true','CoordinateSystem','Geodetic');
    adsbMessage = adsbTx(truePose.Position, truePose.Velocity);

    % Messages can be received within range of each surveillance station
    if isWithinRange(radarsitesLLA, truePose.Position, adsbRange)
        adsbTracks = adsbRx(adsbMessage, time);
    end

    % Update radar tracker at end of scan
    if all([configs.IsValidTime]) && (~isempty(detBuffer) || isLocked(radarGNN))
        radarTracks = radarGNN(detBuffer,time);
        detBuffer = {};
    end
 
    % Fuse tracks
    if isLocked(fuser) || ~isempty([adsbTracks;radarTracks])
        fusedTracks = fuser([adsbTracks;radarTracks],time);
    end

    % plot tracks only once every radar scan
    if all([configs.IsValidTime])

        labels = [repmat(adsblabel,1,numel(adsbTracks)),...
            repmat(radarlabel,1,numel(radarTracks)),...
            repmat(fusedlabel,1,numel(fusedTracks))];
        
        colors = [repmat(adsbclr, numel(adsbTracks), 1) ;...
            repmat(radarclr, numel(radarTracks), 1);...
            repmat(fusedclr, numel(fusedTracks), 1)];
        
        plotTrack(viewer,[adsbTracks; radarTracks; fusedTracks],'ECEF',...
            'LabelStyle','Custom',"CustomLabel",labels,'Color',colors,...
            'LineWidth',3);
    end

    % Record the estimated airplane data for metrics
    fusedTrackLog = [fusedTrackLog, {fusedTracks}]; %#ok<AGROW>
    radarTrackLog = [radarTrackLog, {radarTracks}]; %#ok<AGROW>
    adsbTrackLog = [adsbTrackLog, {adsbTracks}]; %#ok<AGROW>

    % Move camera and take snapshots
    images = moveCamera(viewer,airplane,time,snapTimes,images);

end

figure
imshow(images{1});

В начале сценария самолет далек от самой южной станции наблюдения, и никакое сообщение ADS-B не передается. В результате самолет прослежен радаром только. Обратите внимание на то, что обнаружения ARSR относительно неточны в высоте, которая обычно приемлема как авиадиспетчер отдельные самолеты горизонтально, а не вертикально. Наименее точные обнаружения сгенерированы радарными сайтами, расположенными в более длинных областях значений. Эти обнаружения, тем не менее, сопоставлены к дорожке. На рисунке белая линия представляет истинную траекторию, и желтая линия представляет траекторию, оцененную радарным средством отслеживания. Первый матч рейса далек от любого радара, и соответствующие обнаружения имеют высокий шум измерения. Кроме того, постоянная скоростная модель движения не моделирует движение хорошо во время начального поворота рейса после взлетания. Радарная дорожка передается термофиксатору, который выводит сплавленную дорожку, отображенную оранжевым, после тесно радарной дорожки. Сплавленная дорожка отличается от радарной дорожки, потому что термофиксатор добавляет свой собственный шум процесса в состояние дорожки.

figure
imshow(images{2});

В вышеупомянутом снимке состояния самолет в коммуникационной области значений ADS-B, и устанавливается новая дорожка ADS-B. Термофиксатор обработал этот новый трек и улучшил точность сплавленной дорожки.

figure
imshow(images{3});

В вышеупомянутом снимке состояния самолет вводит конус тишины. Радарное средство отслеживания удаляет дорожку после нескольких обновлений без любых новых обнаружений. В этой точке термофиксатор полагается только на дорожку ADS-B, чтобы оценить положение самолета.

figure
imshow(images{4});

Когда самолет вводит область, покрытую вторыми и третьими радарными станциями наблюдения, новая радарная дорожка устанавливается. Обнаружения от двух радарных станций сплавлены радарным средством отслеживания, и термофиксатор дорожки плавит новую радарную дорожку с дорожкой ADS-B.

Анализ результатов

Вы сравниваете записанные данные о дорожке для радара и ADS-B. Истинное положение и скорость доступны от geoTrajectory. Вы используете метрику OSPA, чтобы сравнить качество отслеживания от радара только, ADS-B только, и от радара, сплавленного с ADS-B. Вы используете настройки по умолчанию trackOSPAMetric объект сравнить NEES (нормированная ошибка расчета придала квадратную форму) для положения дорожки и присвоения обнаружения.

ospa = trackOSPAMetric;
radarospa = zeros(ceil(numel(radarTrackLog)/12),1);

% Compute radar tracker OSPA
for i=1:12:numel(radarTrackLog)
    tracks = radarTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    radarospa(i) = ospa(tracks, truths);
end

% Compute ADS-B OSPA
adsbospa = zeros(numel(adsbTrackLog),1);
reset(ospa);
for i=1:numel(adsbTrackLog)
    tracks = adsbTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    adsbospa(i) = ospa(tracks, truths);
end

% Compute fused OSPA
fusedospa = zeros(numel(fusedTrackLog),1);
reset(ospa);
for i=1:numel(fusedTrackLog)
    tracks = fusedTrackLog{i};
    [truths.Position, ~,truths.Velocity] = lookupPose(flightRoute, i-1,'ECEF');
    fusedospa(i) = ospa(tracks, truths);
end

% Plot OSPA
figure
hold on
plot((0:12:(numel(radarTrackLog)-1))/60, radarospa(1:12:end), "Color",radarclr);
plot((0:(numel(adsbTrackLog)-1))/60,adsbospa, "Color", adsbclr, 'LineWidth',2);
plot((0:(numel(fusedTrackLog)-1))/60,fusedospa, "Color", fusedclr);
l=legend('Radar','ADS-B','Radar + ADS-B');
l.Color = [0.1 0.1 0.1];
l.TextColor = [ 1 1 1];
xlabel('Time (min)')
ylabel('OSPA')
ax = gca;
grid on;
box on;
ax.Color = [0.1 0.1 0.1];
ax.GridColor = [1 1 1];

Метрика OSPA показывает улучшения, полученные из плавления дорожек ADS-B с радарными дорожками. От времени симуляции 19 min к 25 min радар только OSPA высок, потому что самолет пролетает над мертвой точкой радарной сети. Доступность ADS-B в этой области значительно увеличивает эффективность отслеживания, как обозначено сплавленным OSPA. Кроме того, метрика показывает два peaks вначале, которые могут быть приписаны низкой производительности фильтра постоянной скорости во время начальных поворотов траектории и недоступности ADS-B. Около времени 40 min ADS-B только OSPA ухудшается потерей доступности ADS-B в области. В позже segements симуляции, и радар и ADS-B доступны. Радар только OSPA в целом хуже, чем ADS-B только. Это вызвано тем, что радары имеют плохую вертикальную точность по сравнению с GPS.

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

В этом примере вы изучили, как создать сосредоточенный Землей сценарий и задать траектории с помощью геодезических координат. Вы также изучили, как к модели Air Route Surveillance Radars и генерируют синтетические обнаружения. Вы кормите этими обнаружениями мультиобъектное средство отслеживания и оцениваете положение, скорость и заголовок плоскости. Производительность отслеживания улучшается путем добавления и плавления информации о ADS-B. Вы смоделировали отчеты ADS-B и интегрировали их к решению для отслеживания. В этом примере был смоделирован только один рейс. Преимущество ADS-B может быть далее проявлено при моделировании нескольких рейсов в сценариях, где безопасные разделительные расстояния должны быть обеспечены воздушными диспетчерами.

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

initfilter задает расширенный Фильтр Калмана, используемый trackerGNN. Движение самолета хорошо аппроксимировано постоянной скоростной моделью движения. Поэтому используйте относительно маленький шум процесса, чтобы выделить больше веса динамике по сравнению с измерениями, которые, как ожидают, будут довольно шумными в больших расстояниях.

function filter = initfilter(detection)
filter = initcvekf(detection);
filter.StateCovariance = 100*filter.StateCovariance; % initcvekf uses measurement noise as the default
filter.ProcessNoise = 0.1;
end

postProcessDetection постобрабатывает обнаружения двумя операциями:

function detsout = postProcessDetection(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);

velCovTransversal = 100^2;
for i=1:numel(detsout)
    velCov = detsout{i}.MeasurementNoise(4:6,4:6);
    [P,D] = eig(velCov);
    [m,ind] = min(diag(D));
    D = velCovTransversal * eye(3);
    D(ind,ind) = m;
    correctedCov = P*D*P';
    detsout{i}.MeasurementNoise(4:6,4:6) = correctedCov;
end
end

isWithinRange возвращает true, если положение самолета в области значений приемника ADS-B какой-либо из станций наблюдения.

function flag = isWithinRange(stationsLLA, pos, range)
cartDistance = fusion.internal.frames.lla2ecef(stationsLLA) - fusion.internal.frames.lla2ecef(pos);
flag = any(vecnorm(cartDistance,2,2) < range);
end

moveCamera задает новые положения камеры, чтобы следовать за самолетом и взять снимки состояния.

function images = moveCamera(viewer, airplane, time, snapTimes,images)

if mod(time,12) == 0

    pos = airplane.Position;
    if time == 0
        campos(viewer, pos(1),pos(2), 2e4);
    elseif time == 1100 % Zoom out
        campos(viewer,pos(1),pos(2), 1e5);
    elseif time == 2000 % Zoom in
        campos(viewer,pos(1),pos(2), 2.6e4);
    else % Keep previous height but follow airplane
        campos(viewer,pos(1),pos(2));
    end
end

% Snapshots
if any(time == snapTimes)
    drawnow
    img = snapshot(viewer);
    images = [ images, {img}];
end
end

Ссылка

[1] FAR §91.227 Автоматический Зависимый, Широковещательно переданный наблюдением (ADS-B) требования рабочей характеристики