Многоплатформенный радарный Fusion обнаружения

Этот пример показывает, как плавить радарные обнаружения от многоплатформенной радарной сети. Сеть включает два бортовых и наземные радарные платформы дальние. Смотрите пример Multiplatform Radar Detection Generation для деталей. Центральное средство отслеживания обрабатывает обнаружения со всех платформ в фиксированном интервале обновления. Это позволяет вам оценить производительность сети против различных целевых типов и маневров, а также конфигураций платформы и местоположений.

Загрузите моделируемую истину данных и земли датчика

Файл MultiplatformRadarDetectionGeneration.mat содержит три переменные.

  • truthLog является struct, который содержит наземную истину кинематическая информация всех семи целей и платформ на каждом шаге симуляции (сохраненный в поле Truth). Поле Configurations содержит ориентацию, положение луча и кадр координаты датчика для всех четырех радарных датчиков на каждом шаге симуляции.

  • dataLog является struct, который содержит обнаружения, о которых сообщает каждый датчик наряду с неуверенностью и сигналом к шумовому отношению (ОСШ). Времена обнаружения соответствуют временам симуляции по крайней мере с одним обнаружением, о котором сообщают.

  • scene содержит trackingScenario, который используется в целях визуализации.

path = fullfile(matlabroot,'examples','fusion','main');
addpath(path)
load('MultiplatformRadarDetectionGeneration.mat');

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

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

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

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

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

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

Отследите цели путем плавления обнаружений в центральном средстве отслеживания

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

trackUpdateRate = 1/2;   % Update the tracker every 2 seconds
stop = scene.StopTime;
trackUpdateTimes = (1/trackUpdateRate:1/trackUpdateRate:stop)';

% Values used for tracker update logic.
allTracks = [];

% Create a display to show the true, measured, and tracked positions of the
% detected targets and platforms.
theaterDisplay = helperMultiPlatFusionDisplay(scene, 'XLim', [-45 45], ...
    'YLim', [-45 45], 'PlotAssignmentMetrics', true);

detectableTrackIDs = uint32([]);

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

for i=1:numel(trackUpdateTimes)
    time = trackUpdateTimes(i);

    % Buffer all detections since last track update
    isDetReady = (dataLog.DetectionTime>time-1/trackUpdateRate) & (dataLog.DetectionTime<=time);
    detBuffer = dataLog.Detections(isDetReady);

    % Get all beam positions since last track update
    configSelect = (truthLog.SimulationTime>time-1/trackUpdateRate) & (truthLog.SimulationTime<=time);
    configsBuffer = truthLog.Configurations(configSelect,:);

    lastDetectionTime = max(dataLog.DetectionTime(isDetReady));
    if isempty(lastDetectionTime)
        lastDetectionTime = time; % There was no detections since last update
    end

    % Collect list of tracks which fell within at least one radar's field
    % of view since last tracker update
    if isLocked(tracker)
        predictedtracks = predictTracksToTime(tracker, 'all', lastDetectionTime);
        detectableTrackIDs = detectableTracks(allTracks,predictedtracks,configsBuffer);
    end

    % Update tracker, only run track logic on tracks that fell within at
    % least one radar's field of view since last tracker update
    [confirmedTracks, ~, allTracks] = tracker(detBuffer,lastDetectionTime,detectableTrackIDs);

    % Analyze and retrieve the current track-to-truth assignment metrics
    truths = truthLog.Truth(find(configSelect,1,'last'),:);

    [~,~] = tam(confirmedTracks, truths);

    % 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 and tracks
    theaterDisplay(detBuffer,configsBuffer(end,:),confirmedTracks,assignmentTable,truths);

end

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

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

  9x6 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      23               7              20                0                  1                  5        
    60 sec      32             NaN               7                0                  1                  7        
    60 sec      41               3              10                0                  0                  0        

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

maxCondNum = 300;
figure;
helperPlotLongRangeCorrection(maxCondNum)

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

Моделируйте с исправлением ковариации дальним

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

correctedDdets = longRangeCorrection(dataLog.Detections,maxCondNum);
[confirmedTracks,correctedAssignmentTable,ctheaterDisplay] = ...
    runMultiPlatFusionSim(scene,tracker,correctedDdets,dataLog,truthLog);

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

  7x6 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      39              3               10                0                  0                  0        

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

ctheaterDisplay(dataLog.Detections,[],confirmedTracks,correctedAssignmentTable,truths)
axes(ctheaterDisplay.TheaterPlot.Parent)
legend('off')
xlim([-1 5]); ylim([-41 -36]); zlim([-5 0]);
view([-90 90])
axis square
title('Jet Executing Horizontal Turn')

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

view([-60 25])

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

xlim([-25 -9]); ylim([-31 -19]); zlim([-9 -2]);
view([-45 10])
title('Crossing Airliners')

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

xlim([-10 10]); ylim([-0 10]); zlim([-12 -5]);
view([-15 10])
title('Airborne Radar Platforms')
rmpath(path)

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

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

Этот пример показывает, как обработать обнаружения, собранные через несколько бортовых и наземных радарных платформ в центральном средстве отслеживания. В этом примере вы изучили, как шум измерения в больших расстояниях точно не моделируется Распределением Гаусса. Вогнутость эллипса с 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, и структура config, задающая ориентацию датчика относительно координатного кадра дорожки, возвращена объектом 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