Отследите цели точки в плотной помехе Используя средство отслеживания GM-PHD

Этот пример показывает вам, как отследить цели точек в плотной помехе с помощью Гауссовой плотности гипотезы вероятности смеси (GM-PHD) средство отслеживания с помощью постоянной скоростной модели.

Сценарий Setup

Сценарий, используемый в этом примере, создается с помощью trackingScenario. Сценарий состоит из пяти целей точки, перемещающихся в постоянную скорость. Цели перемещаются в поле зрения статического 2D радарного датчика. Вы используете monostaticRadarSensor симулировать 2D датчик и монтировать его на статической платформе. Вы используете FalseAlarmRate свойство датчика управлять плотностью помехи. Значение FalseAlarmRate свойство представляет вероятность генерации ложного предупреждения в одной ячейке разрешения датчика. На основе ложного сигнального уровня 10-3и разрешение датчика, заданного в этом примере, существует приблизительно 48 ложных предупреждений, сгенерированных на шаг.

% Reproducible target locations
rng(2022);

% Create a trackingScenario object
scene = trackingScenario('UpdateRate',10,'StopTime',10);
numTgts = 5;

% Initialize position and velocity of each target
x = 100*(2*rand(numTgts,1) - 1);
y = 100*(2*rand(numTgts,1) - 1);
z = zeros(numTgts,1);
vx = 5*randn(numTgts,1);
vy = 5*randn(numTgts,1);
vz = zeros(numTgts,1);

% Add platforms to scenario with given positions and velocities.
for i = 1:numTgts
    thisTgt = platform(scene);
    thisTgt.Trajectory.Position = [x(i) y(i) z(i)];
    thisTgt.Trajectory.Velocity = [vx(i) vy(i) vz(i)];
end

% Add a detecting platform to the scene.
detectingPlatform = platform(scene);
detectingPlatform.Trajectory.Position = [-200 0 0];

% Simulate 2-D radar using a monostaticRadarSensor.
radar = monostaticRadarSensor(1,...
    'UpdateRate',scene.UpdateRate,...
    'DetectionProbability',0.9,...% Pd
    'FalseAlarmRate',1e-3,...% Pfa
    'FieldOfView',[120 1],...
    'ScanMode','No Scanning',...
    'DetectionCoordinates','scenario',...% Report in scenario frame
    'HasINS',true,...                    % Enable INS to report in scenario frame
    'ReferenceRange',300,...             % Short-range
    'ReferenceRCS',10,...                % Default RCS of platforms
    'MaxUnambiguousRange',400,...        % Short-range
    'RangeResolution',1,...
    'AzimuthResolution',1,...
    'HasFalseAlarms',true);              % Reports false alarms

% Mount the radar on the detecting platform.
detectingPlatform.Sensors = radar;

Настройте средство отслеживания и метрики

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

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

SensorIndex из настройки собирается в 1 совпадать с тем из симулированного датчика. Когда датчик является датчиком точечного объекта, что выходные параметры самое большее одно обнаружение на объект на сканирование, вы устанавливаете MaxNumDetsPerObject свойство настройки к 1.

SensorTransformFcn, SensorTransformParameters, и SensorLimits вместе позвольте вам задавать область, в которой дорожки могут быть обнаружены датчиком. SensorTransformFcn задает преобразование состояния дорожки (xtrack) в промежуточное пространство, использованное датчиком (xsensor) задавать обнаружительную способность дорожки. Полное вычисление, чтобы вычислить вероятность обнаружения показывают ниже:

xsensor=SensorTransformFcn(xtrack,SensorTransformParameters)

Pd=

{Configuration.DetectionProbabilitySensorLimits(:,1)xsensorSensorLimits(:,2)Configuration.MinDetectionProbabilityotherwise

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

Заметьте что подпись SensorTransformFcn похоже на типичную модель измерения. Поэтому можно использовать функции как cvmeas, cameas как SensorTransformFcn. В этом примере вы принимаете, что все дорожки обнаруживаемы. Поэтому SensorTransformFcn задан как @(x,params)x и SensorLimits заданы как [-inf inf] для всех состояний.

% Transform function and limits
sensorTransformFcn = @(x,params)x;
sensorTransformParameters = struct;
sensorLimits = [-inf inf].*ones(6,1);

% Detection probability for sensor configuration
Pd = radar.DetectionProbability;

config = trackingSensorConfiguration('SensorIndex',1,...
    'IsValidTime',true,...% Update every step
    'MaxNumDetsPerObject',1,...
    'SensorTransformFcn',sensorTransformFcn,...
    'SensorTransformParameters',sensorTransformParameters,...
    'SensorLimits',sensorLimits,...
    'DetectionProbability',Pd);

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

% Sensor Parameters to calculate Volume of the resolution bin
Rm = radar.MaxUnambiguousRange/2; % mean -range
dTheta = radar.ElevationResolution;
% Bias fractions reduce the "effective" resolution size
dPhi = radar.AzimuthBiasFraction*radar.AzimuthResolution;
dR = radar.RangeBiasFraction*radar.RangeResolution;

% Cell volume
VCell = 2*((Rm+dR)^3 - Rm^3)/3*(sind(dTheta/2))*deg2rad(dPhi);

% False alarm rate 
Pfa = radar.FalseAlarmRate;

% Define clutter density
config.ClutterDensity = Pfa/VCell;

Вы также задаете FilterInitializationFcn задавать тип фильтра и распределение компонентов в фильтре, инициализированном этим датчиком. В этом примере вы устанавливаете FilterInitializationFcn к initcvgmphd, который создает постоянную скорость фильтр GM-PHD и добавляет 1 компонент на обнаружение низкой вероятности от средства отслеживания. initcvgmphd не добавляет никакой компонент, когда названо без обнаружения. Это означает, что под этой настройкой, компоненты рождения только добавляются к фильтру, когда обнаружения падают за пределами AssignmentThreshold из мультицелевой программы для работы с файлами. См. AssignmentThreshold свойство of trackerPHD для получения дополнительной информации.

config.FilterInitializationFcn = @initcvgmphd;

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

% Use number of radar cells to compute number of false alarms per unit time.
NCells = radar.MaxUnambiguousRange/radar.RangeResolution*radar.FieldOfView(1)/radar.AzimuthResolution;
numFalse = Pfa*NCells;

% Choose an initial weight of 0.05 for each new component. As number of new
% targets is not known, new number of components are simply used as number
% of false alarms. 0.1 is the update rate.
birthRate = 0.05*(numFalse)/0.1;

% Create a tracker with the defined sensor configuration.
tracker = trackerPHD('BirthRate',birthRate,...
    'SensorConfigurations', config);

Метрики

Чтобы оценить производительность средства отслеживания, вы также настраиваете метрику для оценки результатов деятельности. В этом примере вы используете Обобщенное Оптимальное Присвоение Подшаблона (GOSPA) метрика. Метрика GOSPA стремится оценивать производительность средства отслеживания путем присвоения его одна величина затрат. Чем лучше производительность отслеживания, тем ниже GOSPA стоится. Значение нуля представляет совершенное отслеживание.

% Create gospa metric object.
gospa = trackGOSPAMetric;

% Initialize variable for storing metric at each step.
gospaMetric = zeros(0,1);
loc = zeros(0,1);
mt = zeros(0,1);
ft = zeros(0,1);

Запустите симуляцию

Затем вы совершенствуете сценарий, собираете обнаружения из сценария и запускаете средство отслеживания PHD на симулированных обнаружениях.

% Create a display
display = helperClutterTrackingDisplay(scene);

count = 1; % Counter for storing metric data
rng(2018); % Reproducible run

while advance(scene)
    % Current time
    time = scene.SimulationTime;
    
    % Current detections
    detections = detect(scene);
    
    % Tracks
    tracks = tracker(detections, time);
    
    % Update display
    display(scene, detections, tracks);
    
    % Compute GOSPA. getTruth function is defined below.
    truths = getTruth(scene);
    [~,gospaMetric(count),~,loc(count),mt(count),ft(count)] = gospa(tracks, truths);
    count = count + 1;
end

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

figure('Units','normalized','Position',[0.1 0.1 0.8 0.8])
subplot(2,2,[1 2]); plot(gospaMetric,'LineWidth',2); title('Total GOSPA');
ylabel('GOSPA Metric'); xlabel('Time step'); grid on;
subplot(2,2,3); plot(mt,'LineWidth',2); title('Missed Target GOSPA');
ylabel('Missed Target Metric'); xlabel('Time step'); grid on;
subplot(2,2,4); plot(ft,'LineWidth',2); title('False Track GOSPA');
ylabel('False Target Metric'); xlabel('Time step'); grid on;

Анализируйте производительность

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

В этом разделе вы запускаете различную реализацию сценария и средства отслеживания с различными ложными сигнальными уровнями и вычисляете средний GOSPA системы для каждой реализации. Процесс выполнения сценария и вычисления среднего GOSPA для системы перенесен в функции помощника helperRunMonteCarloAnalysis. Чтобы ускорить симуляции Монте-Карло, можно сгенерировать код для monostaticRadarSensor модель, а также trackerPHD использование Тулбокса MATLAB® Coder™. Процесс для генерации кода перенесен в функции помощника helperGenerateCode. Чтобы сгенерировать код для алгоритма, вы собираете код в автономную функцию. Эту функцию называют clutterSimTracker_kernel в этом примере. clutterSimTracker_kernel функция записана, чтобы поддержать четыре ложных сигнальных уровня. Чтобы регенерировать код с различными ложными сигнальными уровнями, можно использовать следующую команду.

Pfas = 10.^(-4,-3,4); % Choose your false alarm settings
helperGenerateCode(scene, Pfas); % Generate code

Чтобы регенерировать код для более ложных сигнальных настроек, можно изменить функции, чтобы включать более персистентные переменные. Для получения дополнительной информации о том, как сгенерировать код для Системы objects™, обратитесь к, Как Сгенерировать код С для примера Средства отслеживания.

helperRunMonteCarloAnalysis функционируйте использует parfor выполнить каждый запущенный Монте-Карло. Если у вас есть лицензия Parallel Computing Toolbox™, запуски Монте-Карло могут быть распределены на нескольких рабочих, чтобы далее ускорить симуляции.

% Turn on flag to run Monte Carlo Simulations
Pfas = 10.^linspace(-4,-3,4);
runMonteCarlo = false;
if runMonteCarlo
    numRunsPerPfa = 50; %#ok<UNRCH> 
    gospaMCs = zeros(4,numRunsPerPfa,4);
    for i = 1:4
        gospaMCs(i,:,:) = helperRunMonteCarloAnalysis(scene, Pfas, Pfas(i), numRunsPerPfa);
    end
    save clutterTrackingMCRuns.mat gospaMCs;
else 
    % Reload results from last run
    load('clutterTrackingMCRuns.mat','gospaMCs');
end

% Plot results
figure('Units','normalized','Position',[0.1 0.1 0.8 0.8])
subplot(2,2,[1 2]); plot(gospaMCs(:,:,1)','LineWidth',2); title('GOSPA');
legend(strcat('Pfa = ',string(Pfas)),'Orientation','horizontal','Location','NorthOutside');
ylabel('Average GOSPA Metric'); xlabel('Monte Carlo Run'); grid on;
subplot(2,2,3); plot(gospaMCs(:,:,3)','LineWidth',2); title('Missed Target GOSPA');
ylabel('Average Missed Target Metric'); xlabel('Monte Carlo Run'); grid on;
subplot(2,2,4); plot(gospaMCs(:,:,4)','LineWidth',2); title('False Track GOSPA');
ylabel('Average False Track Metric'); xlabel('Monte Carlo Run'); grid on;

Графики выше производительности показа средства отслеживания в этом сценарии путем выполнения 50 реализации Монте-Карло для каждого ложного сигнального уровня. Как ложные сигнальные повышения ставки, увеличивается вероятность генерации ложной дорожки. Эта вероятность является четной выше около датчика, где плотность ячеек разрешения намного выше. Когда ложные предупреждения близко расположены и часто появляются в этой области, они могут быть неправильно классифицированы как ложная дорожка низкой скорости. Это поведение средства отслеживания может наблюдаться в усредненном "Ложном Компоненте Дорожки" метрики GOSPA на запущенный сценарий. Заметьте, что как ложные сигнальные повышения ставки, количество peaks в графике также увеличивается. Это также вызывает увеличение общей метрики GOSPA. "Пропущенный Целевой Компонент" является нулем для всех кроме одного запуска. Этот тип события вызывается несколькими промахами цели датчиком.

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

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

Служебные функции

helperRunMonteCarloAnalysis

function gospaMC = helperRunMonteCarloAnalysis(scene, Pfas, Pfai, numRuns)

clear clutterSimTracker_kernel;

% Compute which tracker to run based on provided Pfa
settingToRun = find(Pfas == Pfai,1,'first');

% Initialize ospa for each Monte Carlo run
gospaMC = zeros(numRuns,4);

% Use parfor for parallel simulations
parfor i = 1:numRuns
    rng(i, 'twister');
    
    % Restart scenario before each run
    restart(scene);
    
    % metrics vs time in this scenario
    gospaMetric = zeros(0,4);
    
    % Counter
    count = 0;
    
    % Advance scene and run
    while advance(scene)
        count = count + 1;
        
        % Use radar sensor using targetPoses function
        own = scene.Platforms{end};
        tgtPoses = targetPoses(own,'rotmat');
        insPose = pose(own);
        poses = platformPoses(scene,'rotmat');
        truths = poses(1:end-1);
        time = scene.SimulationTime;
        
        % Reset tracker at the end of step so at the next call tracker is
        % reset to time 0.
        systemToReset = settingToRun*(abs(time - scene.StopTime) < 1e-3);
        
        % Store ospa metric. Tracks and detections can be outputted to run
        % scenario in code generation without needing Monte Carlo runs.
        [~,~,gospaMetric(count,:)] = clutterSimTracker_kernel(Pfas, tgtPoses, insPose, truths, time, settingToRun, systemToReset);
    end
    
    % Compute average of GOSPA in this run
    % Start from time step 20 to allow track establishment.
    gospaMC(i,:) = mean(gospaMetric(20:end,:));
end

end

helperGenerateCode

function helperGenerateCode(scene,Pfas) %#ok<DEFNU>

% Generate sample pose using platformPoses
poses = platformPoses(scene,'rotmat');

% Generate sample inputs for the kernel function
%
% Pfas must be compile-time constant and they cannot change with time
Pfas = coder.Constant(Pfas);

% Target poses as a variable size arrray of max 20 elements
tgtPoses = coder.typeof(poses(1),[20 1],[1 0]);

% Scalar ins pose
insPose = pose(scene.Platforms{1},'true');

% Truth information with maximum 20 targets
truths = coder.typeof(poses(1),[20 1],[1 0]);

% Time
time = 0;

% Which false-alarm setting to run and which to reset. Reset is necessary
% after each run to reinitialize the tracker
systemToRun = 1;
systemToReset = 1;

inputs = {Pfas, tgtPoses, insPose, truths, time, systemToRun, systemToReset}; %#ok<NASGU>

% Same name as MATLAB file allows to shadow the MATLAB function when
% MEX file is available and execute code in MEX automatically.
codegen clutterSimTracker_kernel -args inputs -o clutterSimTracker_kernel;

end

getTruth

function truths = getTruth(scenario)
platPoses = platformPoses(scenario); % True Information
truths = platPoses(1:end-1); % Last object is ownship
end

clutterSimTracker_kernel

Эта функция задана во внешнем файле с именем clutterSimTracker_kernel, доступный в той же рабочей папке как этот скрипт.

function [detections, tracks, ospaMetric] = clutterSimTracker_kernel(Pfas, tgtPoses, insPose, truths, time, systemToRun, systemToReset)

assert(numel(Pfas) == 4,'Only 4 false alarm settings supported. Rewrite more persistent variables to add more settings');

persistent tracker1 tracker2 tracker3 tracker4 ...
    radar1 radar2 radar3 radar4 ....
    reset1 reset2 reset3 reset4 ....
    gospa

if isempty(gospa) || isempty(reset1) || isempty(reset2) || isempty(reset3) || isempty(reset4)
    gospa = trackGOSPAMetric('CutoffDistance',50);
end

if isempty(tracker1) || isempty(radar1) || isempty(reset1)
    tracker1 = setupTracker(Pfas(1));
    radar1 = setupRadar(Pfas(1));
    reset1 = zeros(1,1);
end

if isempty(tracker2) || isempty(radar2) || isempty(reset2)
    tracker2 = setupTracker(Pfas(2));
    radar2 = setupRadar(Pfas(2));
    reset2 = zeros(1,1);
end

if isempty(tracker3) || isempty(radar3) || isempty(reset3)
    tracker3 = setupTracker(Pfas(3));
    radar3 = setupRadar(Pfas(3));
    reset3 = zeros(1,1);
end

if isempty(tracker4) || isempty(radar4) || isempty(reset4)
    tracker4 = setupTracker(Pfas(4));
    radar4 = setupRadar(Pfas(4));
    reset4 = zeros(1,1);
end

switch systemToRun
    case 1
        detections = radar1(tgtPoses, insPose, time);
        tracks = tracker1(detections, time);
    case 2
        detections = radar2(tgtPoses, insPose, time);
        tracks = tracker2(detections, time);
    case 3
        detections = radar3(tgtPoses, insPose, time);
        tracks = tracker3(detections, time);
    case 4
        detections = radar4(tgtPoses, insPose, time);
        tracks = tracker4(detections, time);
    otherwise
        error('Idx out of bounds');
end

[~,gp,~,loc,mT,fT] = gospa(tracks, truths);

ospaMetric = [gp loc mT fT];

switch systemToReset
    case 1
        reset1 = zeros(0,1);
    case 2
        reset2 = zeros(0,1);
    case 3
        reset3 = zeros(0,1);
    case 4
        reset4 = zeros(0,1);
end

end

function tracker = setupTracker(Pfa)
Pd = 0.9;
VCell = 0.0609; % Inlined value of cell volume;
Kc = Pfa/VCell; % Clutter density

% Birth rate
birthRate = 0.05*(5 + Pfa*48000)/0.1;

% Transform function and limits
sensorTransformFcn = @(x,params)x;
sensorTransformParameters = struct;
sensorLimits = bsxfun(@times,[-inf inf],ones(6,1));

% Assemble information into a trackingSensorConfiguration
config = trackingSensorConfiguration('SensorIndex',1,...
    'IsValidTime',true,...% Update every step
    'DetectionProbability',Pd,...% Detection Probability of detectable states
    'ClutterDensity',Kc,...% Clutter Density
    'SensorTransformFcn',sensorTransformFcn,...
    'SensorTransformParameters',sensorTransformParameters,...
    'SensorLimits',sensorLimits,...
    'MaxNumDetsPerObject',1,...
    'FilterInitializationFcn',@initcvgmphd);

tracker = trackerPHD('SensorConfigurations', config,...
    'BirthRate',birthRate);

end

function radar = setupRadar(Pfa)
Pd = 0.9;
radar = monostaticRadarSensor(1,...
    'UpdateRate',10,...
    'DetectionProbability',Pd,...
    'FalseAlarmRate',Pfa,...
    'FieldOfView',[120 1],...
    'ScanMode','No Scanning',...
    'DetectionCoordinates','Scenario',...% Report in scenario frame
    'HasINS',true,...% Enable INS to report in scenario frame
    'ReferenceRange',300,...% Short-range
    'ReferenceRCS',10,...
    'MaxUnambiguousRange',400,...% Short-range
    'RangeResolution',1,...
    'AzimuthResolution',1,...
    'HasFalseAlarms',true);
end
Для просмотра документации необходимо авторизоваться на сайте