В этом примере показано, как отслеживать целевые точки в плотном нагромождении с помощью отслеживания плотности вероятностей гауссовой смеси (GM-PHD) с использованием модели постоянной скорости.
Сценарий, используемый в этом примере, создается с помощью trackingScenario. Сценарий состоит из пяти точечных целей, движущихся с постоянной скоростью. Цели перемещаются в поле зрения статического 2-D радиолокационного датчика. Вы используете monostaticRadarSensor имитировать датчик 2-D и монтировать его на статической платформе. Вы используете FalseAlarmRate свойство датчика контролировать плотность беспорядка. Значение FalseAlarmRate свойство представляет вероятность генерации ложного аварийного сигнала в одной ячейке разрешения датчика. На основании частоты ложных аварийных сигналов и разрешения датчика, определенного в этом примере, на шаг генерируется приблизительно 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 определяет преобразование состояния дорожки () в промежуточное пространство, используемое датчиком () для определения возможности обнаружения дорожки. Общий расчет для вычисления вероятности обнаружения показан ниже:
SensorTransformParameters)
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 многозначного файлера. См. АssignmentThreshold имущество 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 для системы переносится в вспомогательную функцию helperRunteCarloAnalysis. Для ускорения моделирования Монте-Карло можно создать код для monostaticRadarSensor модель, а также trackerPHD с помощью MATLAB ® Coder™ Toolbox. Процесс создания кода заключен в вспомогательную функцию helperGenerateCode. Чтобы создать код для алгоритма, необходимо собрать код в автономную функцию. Эта функция названаclutterSimTracker_kernel в этом примере. clutterSimTracker_kernel функция записывается для поддержки четырех частот ложных аварийных сигналов. Для регенерации кода с различными частотами ложных аварийных сигналов можно использовать следующую команду.
Pfas = 10.^(-4,-3,4); % Choose your false alarm settings helperGenerateCode(scene, Pfas); % Generate code
Для регенерации кода для дополнительных параметров ложных аварийных сигналов можно изменить функции, включив в них более постоянные переменные. Дополнительные сведения о создании кода для системного objects™ см. в примере «Создание кода C для трекера».
helperRunMonteCarloAnalysis функция использует parfor для выполнения каждого прогона Монте-Карло. При наличии лицензии Parallel Computing Toolbox™ программы Monte Carlo могут быть распределены между несколькими работниками для дальнейшего ускорения моделирования.
% 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 для каждого прогона сценария. Обратите внимание, что по мере увеличения частоты ложных аварийных сигналов количество пиков на графике также увеличивается. Это также вызывает увеличение общей метрики GOSPA. «Пропущенный целевой компонент» равен нулю для всех, кроме одного прогона. Этот тип события вызван несколькими промахами цели датчиком.
В этом примере вы научились настраивать и инициализировать трекер GM-PHD для отслеживания целевых точек для заданной частоты ложных аварийных сигналов. Также вы научились оценивать производительность трекера с помощью метрики GOSPA и связанных с ней компонентов. Кроме того, вы научились выполнять несколько реализаций сценария при различных настройках ложных аварийных сигналов для определения характеристик производительности трекера.
helperRunMonteCarloAnalysisfunction 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
helperGenerateCodefunction 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
getTruthfunction 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