Загрузка записи сценария отслеживания

The MultiplatformRadarDetectionGeneration MAT-файл содержит запись сценария отслеживания, ранее сгенерированную с помощью следующей команды

recording = record(scene,'IncludeSensors',true,'InitialSeed',2018,'RecordingFormat','Recording')

где scene- сценарий отслеживания, созданный в примере Multiplatform Radar Detection Generation.

load('MultiplatformScenarioRecording.mat');

Задайте центральный трекер

Используйте trackerGNN как центральный трекер, который обрабатывает обнаружения, полученные со всех радиолокационных платформ в сценарии.

Трекер использует initFilter поддерживающая функция для инициализации расширенного фильтра Калмана с постоянной скоростью для каждой новой дорожки. initFilter изменяет фильтр, возвращенный initcvekf чтобы соответствовать высоким целевым скоростям. Технологический шум фильтра установлен в 1$\mathit{g}$ ($\mathit{g}=9.8\mathrm{m}/{\mathrm{s}}^2$) обеспечить сопровождение маневрирующих целей в сценарии.

Трекер AssignmentThreshold установлено равным 50, чтобы позволить обнаружениям с большим смещением области значений (из-за эффектов атмосферного преломления в длинных диапазонах обнаружения) быть сопоставленными с треками в трекере. The DeletionThreshold установлено значение 3, чтобы быстро удалить избыточные дорожки.

Включите HasDetectableTrackIDsInput для определения дорожек, которые находятся в поле зрения по меньшей мере одного радара с момента последнего обновления. Логика трека оценивается только на треках, у которых была возможность обнаружения с момента последнего обновления трекера.

tracker = trackerGNN('FilterInitializationFcn', @initFilter, ...
        'AssignmentThreshold', 50, 'DeletionThreshold', 3, ...
        'HasDetectableTrackIDsInput', true);

Отслеживайте цели при помощи Fusing Detections в центральном трекере

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

trackUpdateRate = 0.5;   % Update the tracker every 2 seconds

% Create a display to show the true, measured, and tracked positions of the
% detected targets and platforms.
theaterDisplay = helperMultiPlatFusionDisplay(recording,'PlotAssignmentMetrics', true);

% Construct an object to analyze assignment and error metrics
tam = trackAssignmentMetrics('DistanceFunctionFormat','custom',...
    'AssignmentDistanceFcn',@truthAssignmentDistance,...
    'DivergenceDistanceFcn',@truthAssignmentDistance);

% Initialize buffers
detBuffer = {};
sensorConfigBuffer = [];
allTracks = [];
detectableTrackIDs = uint32([]);
assignmentTable = [];

% Initialize next tracker update time
nextTrackUpdateTime = 2;

while ~isDone(recording)
    % Read the next record of the recording.
    [time, truePoses, covcon, dets, senscon, sensPlatIDs] = read(recording);
    
    % Buffer all detections and sensor configurations until the next tracker update
    detBuffer = [detBuffer ; dets]; %#ok<AGROW>
    sensorConfigBuffer = [sensorConfigBuffer ; senscon']; %#ok<AGROW>
        
    % Follow the trackUpdateRate to update the tracker
    if time >= nextTrackUpdateTime || isDone(recording)
        
        if isempty(detBuffer)
            lastDetectionTime = time;
        else
            lastDetectionTime = detBuffer{end}.Time;
        end
        
        if isLocked(tracker)
            % Collect list of tracks which fell within at least one radar's field
            % of view since the last tracker update
            predictedtracks = predictTracksToTime(tracker, 'all', lastDetectionTime);
            detectableTrackIDs = detectableTracks(allTracks, predictedtracks, sensorConfigBuffer);
        end
        
        % Update tracker. Only run track logic on tracks that fell within at
        % least one radar's field of view since the last tracker update
        [confirmedTracks, ~, allTracks] = tracker(detBuffer, lastDetectionTime, detectableTrackIDs);
        
        % Analyze and retrieve the current track-to-truth assignment metrics
        tam(confirmedTracks, truePoses);
        
        % Store assignment metrics in a table
        currentAssignmentTable = trackMetricsTable(tam);
        rowTimes = seconds(time*ones(size(currentAssignmentTable,1),1));
        assignmentTable = cat(1,assignmentTable,table2timetable(currentAssignmentTable,'RowTimes',rowTimes));
        
        % Update display with detections, coverages, and tracks
        theaterDisplay(detBuffer, covcon, confirmedTracks, assignmentTable, truePoses, sensPlatIDs);
        
        % Empty buffers
        detBuffer = {};
        sensorConfigBuffer = [];
        
        % Update next track update time
        nextTrackUpdateTime = nextTrackUpdateTime + 1/trackUpdateRate;
    end
    
end

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

endTime = assignmentTable.Time(end);
assignmentTable(endTime,{'TrackID','AssignedTruthID','TotalLength','DivergenceCount','RedundancyCount','RedundancyLength'})

ans=9×6 timetable
     Time     TrackID    AssignedTruthID    TotalLength    DivergenceCount    RedundancyCount    RedundancyLength
    ______    _______    _______________    ___________    _______________    _______________    ________________

    60 sec       1               2              27                0                  0                  0        
    60 sec       7               4              27                0                  0                  0        
    60 sec       8               5              26                0                  0                  0        
    60 sec       9               1              22                0                  0                  0        
    60 sec      11             NaN              10                0                  0                  0        
    60 sec      12               6              20                0                  0                  0        
    60 sec      24               7              19                0                  1                  4        
    60 sec      32             NaN               7                0                  1                  7        
    60 sec      41               3              10                0                  0                  0        

Заметьте, что платформы, которые испытывают трудности с обслуживанием путей (платформы 4 и 7), также являются платформами, наиболее удаленными от радаров. Эта низкая эффективность отслеживания объясняется допущением Гауссова распределения для шума измерения. Предположение хорошо работает для целей с короткими областями значений, но при больших областях значений неопределенность измерения отклоняется от Гауссова распределения. Следующий рисунок сравнивает эллипсы ковариации 1-сигма, соответствующие фактическому распределению цели и распределению цели, заданному радарным датчиком. Датчик находится на расстоянии 5 км от цели с угловым разрешением 5 степеней. Фактическая неопределенность измерения имеет вогнутую форму, вытекающую из координатной системы координат обнаружения сферического датчика, в которой радар оценивает положение цели.

maxCondNum = 300;
figure;
helperPlotLongRangeCorrection(maxCondNum)

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

Симулируйте с большой ковариационной коррекцией

Повторите предыдущую симуляцию с помощью longRangeCorrection поддерживающая функция для коррекции сообщенного шума измерения на больших областях значений.

[confirmedTracks,correctedAssignmentTable,ctheaterDisplay] = ...
    runMultiPlatFusionSim(recording,tracker,@longRangeCorrection);

endTime = correctedAssignmentTable.Time(end);
correctedAssignmentTable(endTime,{'TrackID','AssignedTruthID','TotalLength','DivergenceCount','RedundancyCount','RedundancyLength'})

ans=7×6 timetable
     Time     TrackID    AssignedTruthID    TotalLength    DivergenceCount    RedundancyCount    RedundancyLength
    ______    _______    _______________    ___________    _______________    _______________    ________________

    60 sec       1              2               27                0                  0                  0        
    60 sec       7              4               27                0                  0                  0        
    60 sec       8              5               26                0                  0                  0        
    60 sec       9              1               22                0                  0                  0        
    60 sec      11              7               25                0                  0                  0        
    60 sec      12              6               20                0                  0                  0        
    60 sec      38              3               10                0                  0                  0        

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

allDetections = vertcat(recording.RecordedData.Detections);
ctheaterDisplay(allDetections,covcon,confirmedTracks,correctedAssignmentTable,truePoses, sensPlatIDs) 
axes(ctheaterDisplay.TheaterPlot.Parent)
legend('off')
xlim([-1000 5000]); ylim([-41000 -36000]); zlim([-5000 0]);
view([-90 90])
axis square
title('Jet Executing Horizontal Turn')

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

view([-60 25])

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

xlim([-25000 -9000]); ylim([-31000 -19000]); zlim([-9000 -2000]);
view([-45 10])
title('Crossing Airliners')

Переключая точку зрения, чтобы сосредоточиться на пересекающих лайнерах, изображены те же неточные измерения высоты. Обратите внимание, как красные обнаружения центрируются на высоте 8 км, в то время как два лайнера летают на высотах 3 и 4 км соответственно. Использование очень большой ковариации на высоте позволяет трекеру игнорировать ошибочные показания высоты от красных обнаружений и отслеживать высоту с помощью двух других радаров. Наблюдая ковариацию неопределенности треков T07 и T08, можно увидеть, что они дают последовательную оценку платформ P04 и P05, соответственно.

xlim([-10000 10000]); ylim([-0 10000]); zlim([-12000 -5000]);
view([-15 10])
title('Airborne Radar Platforms')

Последний график фокусируется на двух радиолокационных платформах воздушного базирования. Каждая платформа обнаруживается другой платформой, а также наземным радаром. Траектории платформы пересекают друг друга, разделяясь на 1000 м в высоту, и их пути соответствуют основной истине.

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

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

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

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

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 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 allow for some horizontal maneuver
scaleAccel = cast(10,classToUse);
Q = blkdiag(scaleAccel^2, scaleAccel^2, 1);
filter.ProcessNoise = Q;
end

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

function trackIDs = detectableTracks(tracks,predictedtracks,configs)
% Identify which tracks fell within a sensor's field of view

numTrks = size(tracks,1);
[numsteps, numSensors] = size(configs);
allposTrack = zeros(3,numsteps);
isDetectable = false(numTrks,1);
for iTrk = 1:numTrks
    % Interpolate track positions between current position and predicted
    % positions for each simulation step
    posTrack = tracks(iTrk).State(1:2:end);
    posPreditedTrack = predictedtracks(iTrk).State(1:2:end);
    for iPos = 1:3
        allposTrack(iPos,:) = linspace(posTrack(iPos),posPreditedTrack(iPos),numsteps);
    end
    for iSensor = 1:numSensors
        thisConfig = configs(:,iSensor);
        for k = 1:numsteps
            if thisConfig(k).IsValidTime
                pos = trackToSensor(allposTrack(:,k),thisConfig(k));
                % Check if predicted track position is in sensor field of
                % view
                [az,el] = cart2sph(pos(1),pos(2),pos(3));
                az = az*180/pi;
                el = el*180/pi;
                inFov = abs(az)<thisConfig(k).FieldOfView(1)/2 && abs(el) < thisConfig(k).FieldOfView(2)/2;
                if inFov
                    isDetectable(iTrk) = inFov;
                    k = numsteps; %#ok<FXSET>
                    iSensor = numSensors; %#ok<FXSET>
                end
            end
        end
    end

end

trackIDs = [tracks.TrackID]';
trackIDs = trackIDs(isDetectable);
end

trackToSensor Эта функция возвращает положение дорожки в координатной системе координат датчика. Структура дорожки возвращается trackerGNN объект и конфигурационная структура, определяющая ориентацию датчика относительно координатной системы координат дорожки, возвращаются monostaticRadarSensor объект.

function pos = trackToSensor(pos,config)

frames = config.MeasurementParameters;
for m = numel(frames):-1:1
    rotmat = frames(m).Orientation;
    if ~frames(m).IsParentToChild
        rotmat = rotmat';
    end
    offset = frames(m).OriginPosition;
    pos = bsxfun(@minus,pos,offset);
    pos = rotmat*pos;
end
end

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

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

function dets = longRangeCorrection(dets,maxCond)
for m = 1:numel(dets)
    R = dets{m}.MeasurementNoise;
    [Q,D] = eig(R);
    Q = real(Q);
    d = real(diag(D));
    dMax = max(d);
    condNums = dMax./d;
    iFix = condNums>maxCond;
    d(iFix) = dMax/maxCond;
    R = Q*diag(d)*Q';
    dets{m}.MeasurementNoise = R;
end
end