利用GM-PHD跟踪器跟踪密集杂波中的点目标

打开生活的脚本

这个例子向您展示了如何使用一个恒定速度模型的高斯混合概率假设密度(GM-PHD)跟踪器来跟踪密集杂波中的点目标。

设置场景

本示例中使用的场景是使用trackingScenario．该场景由5个以恒定速度移动的点目标组成。目标在静态二维雷达传感器的视场内移动。您使用monostaticRadarSensor模拟二维传感器，并将其安装在一个静态平台上。您使用FalseAlarmRate利用传感器的特性来控制杂波密度。价值FalseAlarmRate属性表示在传感器的一个分辨单元中产生假警报的概率。基于虚警率 $10^{- 3.}$ 根据本例中定义的传感器的分辨率，每步大约产生48个误报。

可重复的目标位置rng (2022);创建一个trackingScenario对象现场= trackingScenario (“UpdateRate”10“StopTime”10);numTgts = 5;初始化每个目标的位置和速度x = 100*(2*rand(numTgts,1) - 1);y = 100*(2*rand(numTgts,1) - 1);z = 0 (numTgts, 1);vx = 5 * randn (numTgts, 1);v = 5 * randn (numTgts, 1);vz = 0 (numTgts, 1);%在给定位置和速度的情况下添加平台。为i = 1:numTgts thisTgt =平台(场景);thisTgt.Trajectory.Position = [x(i) y(i) z(i)];this . trajectory . velocity = [vx(i) vy(i) vz(i)];结束在场景中添加一个检测平台。detectingPlatform =平台(现场);detection platform . trajectory . position = [-200 00];使用单站式传感器模拟二维雷达。雷达= monostaticRadarSensor (1,.．.“UpdateRate”场景。UpdateRate,.．.“DetectionProbability”, 0.9,.．.% Pd“FalseAlarmRate”1 e - 3,.．.% Pfa“FieldOfView”, 120年[1],.．.“ScanMode”,“没有扫描”,.．.“DetectionCoordinates”,“场景”,.．.%在场景框架中报告“之内”,真的,.．.%启用INS在场景框架中报告“ReferenceRange”, 300,.．.%短程“ReferenceRCS”10.．.%平台默认RCS“MaxUnambiguousRange”, 400,.．.%短程“RangeResolution”,1，.．.“AzimuthResolution”,1，.．.“HasFalseAlarms”,真正的);报告误报将雷达安装在探测平台上。detectingPlatform。传感器=雷达;

设置跟踪器和指标

跟踪器

使用GM-PHD点对象跟踪器来跟踪目标。配置PHD跟踪器的第一步是设置传感器的配置。使用trackingSensorConfiguration对象。

的SensorIndex的配置设置为1，以匹配模拟传感器。由于传感器是点目标传感器，每次扫描每个目标最多输出一次检测，因此设置MaxNumDetsPerObject属性的配置为1。

的SensorTransformFcn,SensorTransformParameters,SensorLimits一起允许您定义传感器可以检测轨道的区域。的SensorTransformFcn定义轨道状态的转换( $x_{跟踪}$ )进入传感器使用的中间空间( $x_{传感器}$ )来定义轨迹可检测性。计算检测概率的总体计算如下图所示:

$x_{传感器} ＝ SensorTransformFcn （ x_{跟踪}, SensorTransformParameters ）$

$P_{d} ＝$

$｛ \begin{array}{l} 配置． DetectionProbability & {SensorLimits （：, 1 ） \leq x}_{传感器} \leq SensorLimits （：, 2 ） \\ 配置． MinDetectionProbability & 否则 \end{array}$

为了计算给定状态协方差的不确定状态的检测概率，跟踪器使用类似于无迹卡尔曼滤波的sigma点计算生成状态样本。

注意SensorTransformFcn类似于一个典型的度量模型。因此，可以像这样使用函数cvmeas,cameas作为SensorTransformFcn．在本例中，假设所有轨迹都是可检测的。因此,SensorTransformFcn被定义为@ x (x, params)和SensorLimits的定义为(负无穷到正无穷)所有国家。

变换函数和极限sensorTransformFcn = @ x (x, params);sensorTransformParameters =结构;sensorLimits = [-inf inf].*ones(1,1);%传感器配置的检测概率Pd = radar.DetectionProbability;配置= trackingSensorConfiguration (“SensorIndex”,1，.．.“IsValidTime”,真的,.．.%更新每一步“MaxNumDetsPerObject”,1，.．.“SensorTransformFcn”sensorTransformFcn,.．.“SensorTransformParameters”sensorTransformParameters,.．.“SensorLimits”sensorLimits,.．.“DetectionProbability”, Pd);

的ClutterDensity配置的属性是指测量空间的单位体积的误报率。在本例中，在场景框架中报告检测时，测量空间被定义为笛卡尔坐标。由于传感器在笛卡尔坐标系下的分辨率体积随分辨率的方位角和距离的变化，可以在传感器的平均距离处计算出一个近似值。

%传感器参数来计算分辨率仓的体积Rm = radar.MaxUnambiguousRange / 2;%的意思是范围dTheta = radar.ElevationResolution;%偏置分数减少了“有效”分辨率大小dPhi = radar.AzimuthBiasFraction * radar.AzimuthResolution;= radar.RangeBiasFraction * radar.RangeResolution博士;细胞体积百分比VCell = 2 * ((Rm +博士)^ 3 - Rm ^ 3) / 3 *(信德(dTheta / 2)) *函数(dPhi);误报率%Pfa = radar.FalseAlarmRate;定义杂波密度配置。ClutterDensity = Pfa / VCell;

你还可以定义FilterInitializationFcn指定过滤器的类型和过滤器中组件的分布，由该传感器初始化。在本例中，您设置FilterInitializationFcn来initcvgmphd，它创建一个恒定速度的GM-PHD滤波器，并从跟踪器中为每一个低可能性检测添加一个组件。的initcvgmphd在没有检测的情况下调用时不添加任何组件。这意味着在这种配置下，只有当检测不在时，出生组件才会被添加到过滤器中AssignmentThreshold多目标的纳税人。看到一个ssignmentThreshold的属性trackerPHD为更多的细节。

配置。FilterInitializationFcn = @initcvgmphd;

接下来，使用该配置创建跟踪器trackerPHD系统对象™。在配置跟踪器时，指定出生率属性来定义单位时间内出现在视野中的目标数量。的FilterInitializationFcn与配置一起使用，为每个未分配的检测添加一个组件。在每个时间步骤中，您可以预期组件的数量大致等于错误警报和新目标的数量。跟踪器分配出生率对所有这些分量都是一致的。

%使用雷达单元数计算单位时间内的误报次数。ncell = radar.MaxUnambiguousRange / radar.RangeResolution * radar.FieldOfView (1) / radar.AzimuthResolution;numFalse = Pfa * ncell;为每个新组件选择0.05的初始权重。如新%目标是未知的，新组件的数量简单地用作数字假警报百分比。0.1为更新速率。出生率= 0.05 * (numFalse) / 0.1;使用定义的传感器配置创建跟踪器。追踪= trackerPHD (“出生率”出生率,.．.“SensorConfigurations”、配置);

指标

要评估跟踪器的性能，还需要设置性能评估指标。在本例中，您将使用广义最优子模式分配(GOSPA)指标。GOSPA指标旨在通过指定一个成本值来评估跟踪器的性能。跟踪性能越好，GOSPA成本越低。0代表完美的跟踪。

创建gospa度量对象。gospa = trackGOSPAMetric;%初始化变量，用于在每一步存储度量。gospaMetric = 0 (0, 1);loc = 0 (0, 1);太= 0 (0,1);英国《金融时报》= 0 (0,1);

运行仿真

接下来，推进该场景，从该场景收集检测，并在模拟检测上运行PHD跟踪器。

创建一个显示显示= helperClutterTrackingDisplay(现场);数= 1;%用于存储度量数据的计数器rng (2018);%的运行而推进(场景)%当前时间时间= scene.SimulationTime;%电流检测检测=检测(现场);%的踪迹Tracks =跟踪器(检测，时间);%更新显示显示器(场景,检测、跟踪);%计算GOSPA。getTruth函数定义如下。真理= getTruth(现场);[~，gospaMetric(count)，~，loc(count)，mt(count)，ft(count)] = gospa(tracks, truths);Count = Count + 1;结束

在这个场景中，您可以看到PHD跟踪器跟踪了所有目标。这也可以使用GOSPA指标和相关组件进行定量评估。在下面的图中，注意GOSPA指标在几个步骤之后下降了。由于每个轨道的建立延迟，GOSPA度量的初始值较高。

图(“单位”,“归一化”,“位置”，[0.1 0.1 0.8 0.8])子plot(2,2，[1 2]);情节(gospaMetric“线宽”2);标题(“总GOSPA”）;ylabel (“GOSPA指标”）;包含(“时间步”）;网格在;次要情节(2、2、3);情节(太“线宽”2);标题(“错过了目标GOSPA”）;ylabel (“错过了目标指标”）;包含(“时间步”）;网格在;次要情节(2、2、4);图(见英国《金融时报》“线宽”2);标题(“虚假跟踪GOSPA”）;ylabel (“虚假目标指标”）;包含(“时间步”）;网格在;

分析性能

确定跟踪器性能的一个典型方法是对场景的不同实现运行几个模拟。蒙特卡罗模拟有助于消除随机事件的影响，如假警报的位置、目标脱靶事件和测量中的噪声。

在本节中，您将运行场景的不同实现和具有不同误报率的跟踪器，并计算每个实现的系统的平均GOSPA。运行场景和计算系统平均GOSPA的过程被包装在helper函数中helperRunMonteCarloAnalysis．为了加速蒙特卡罗模拟，您可以生成代码monostaticRadarSensor模型以及trackerPHD使用MATLAB®编码器™工具箱。生成代码的过程被包装在helper函数中helperGenerateCode．要为算法生成代码，您需要将代码组装到一个独立的函数中。这个函数被命名为clutterSimTracker_kernel在这个例子中。的clutterSimTracker_kernel函数支持四种误报率。金宝app要以不同的误报率重新生成代码，可以使用以下命令。

pfa = 10 ^(4、3、4);选择错误警报设置helperGenerateCode(场景,pfa);%生成代码

要重新生成代码以进行更多的错误警报设置，可以修改函数以包含更多持久变量。有关如何为System对象™生成代码的详细信息，请参阅如何为跟踪器生成C代码的例子。

的helperRunMonteCarloAnalysis函数使用一个parfor来执行每一次蒙特卡洛运行。如果您有Parallel Computing Toolbox™许可证，那么蒙特卡洛运行可以分布到多个工作程序，以进一步加速模拟。

%开启标志以运行蒙特卡罗模拟pfa = 10。^ linspace(4、3、4);runMonteCarlo = false;如果runMonteCarlo numRunsPerPfa = 50;% #好< UNRCH >gospaMCs = 0 (4 numRunsPerPfa 4);为i = 1:4 gospaMCs(i，:，:) = helperRunMonteCarloAnalysis(scene, Pfas, Pfas(i)， numRunsPerPfa);结束保存clutterTrackingMCRuns.matgospaMCs;其他的重新加载上次运行的结果负载(“clutterTrackingMCRuns.mat”,“gospaMCs”）;结束%绘制结果图(“单位”,“归一化”,“位置”，[0.1 0.1 0.8 0.8])子plot(2,2，[1 2]);情节(gospaMCs(:: 1)”,“线宽”2);标题(“GOSPA”）;传奇(strcat (“Pfa = '字符串(pfa)),“定位”,“水平”,“位置”,“NorthOutside”）;ylabel (“平均GOSPA指标”）;包含(“蒙特卡罗运行”）;网格在;次要情节(2、2、3);情节(gospaMCs(:: 3)”,“线宽”2);标题(“错过了目标GOSPA”）;ylabel (“平均未完成目标指标”）;包含(“蒙特卡罗运行”）;网格在;次要情节(2、2、4);情节(gospaMCs(:: 4)”,“线宽”2);标题(“虚假跟踪GOSPA”）;ylabel (“平均虚假轨迹指标”）;包含(“蒙特卡罗运行”）;网格在;

上面的图显示了跟踪器在这种情况下的性能，对每个虚警率运行50个蒙特卡罗实现。随着虚警率的增加，产生虚假航迹的概率也增加。这种可能性在传感器附近更高，那里的分辨率单元密度更高。由于该区域的虚警分布较密集，且频繁出现，容易被误分类为低速伪迹。跟踪器的这种行为可以在每个场景运行的GOSPA度量的平均“假跟踪组件”中观察到。注意，随着虚警率的增加，图中峰值的数量也会增加。这也会导致总的GOSPA指标的增加。除了一次运行外，“错过的目标组件”对所有运行都是零。这种类型的事件是由传感器多次未命中目标引起的。

总结

在本例中，您了解了如何配置和初始化GM-PHD跟踪器，以跟踪给定假警报率的点目标。您还了解了如何使用GOSPA指标及其相关组件评估跟踪器的性能。此外，您还学习了如何在不同的假警报设置下运行场景的几个实现，以确定跟踪器的性能特征。

效用函数

`helperRunMonteCarloAnalysis`

函数gospaMC = helperRunMonteCarloAnalysis(scene, Pfas, Pfai, numRuns) clearclutterSimTracker_kernel;%根据提供的Pfa计算要运行的跟踪器setingtorun = find(pas == Pfai,1，)“第一”）;%为每次蒙特卡洛运行初始化ospagospaMC = 0 (numRuns 4);%使用parfor进行并行模拟parfori = 1:numRuns rng(i，“旋风”）;每次运行前重新启动场景重启(现场);%度量与时间gospaMetric = 0 (0, 4);%计数器数= 0;%前进场景和运行而Advance (scene) count = count + 1;使用targetpose功能使用雷达传感器结束自己= scene.Platforms {};tgtPoses = targetPoses(自己的,“rotmat”）;insPose =姿势(自己的);提出了= platformPoses(场景,“rotmat”）;真理=姿势(1:end-1);时间= scene.SimulationTime;在步骤结束时重置跟踪器，以便在下一个调用跟踪器时重置时间0。systemToReset = settingToRun*(abs(time - scene.StopTime) < 1e-3); / /设定时间%存储ospa度量。轨迹和检测可以输出运行%场景的代码生成，而不需要运行蒙特卡罗。[~，~，gospaMetric(count，:)] = clutterSimTracker_kernel(Pfas, tgtPoses, insPose, truths, time, settingToRun, systemToReset);结束%计算这次运行中GOSPA的平均值%从时间步骤20开始，允许建立跟踪。gospaMC(我)=意味着(gospaMetric(20:最终,));结束结束

`helperGenerateCode`

函数helperGenerateCode(场景,pfa)% #好< DEFNU >使用platformpose生成样本姿势提出了= platformPoses(场景,“rotmat”）;%为内核函数生成样本输入％% Pfas必须是编译时常量，并且不能随时间而改变pfa = coder.Constant (pfa);% Target作为一个不超过20个元素的可变大小数组tgtpose = code .typeof(pose (1)，[20 1]，[1 0]);%标量ins位姿insPose =姿势(场景。平台{1},“真正的”）;%真实信息，最多20个目标truth = code .typeof(pose (1)，[20 1]，[1 0]);%的时间时间= 0;要运行哪个假警报设置，哪个要复位。重置是必要的%在每次运行后重新初始化跟踪器systemToRun = 1;systemToReset = 1;input = {Pfas, tgtpose, insPose, truths, time, systemToRun, systemToReset};% #好< NASGU >与MATLAB同名的文件允许阴影MATLAB函数时% MEX文件可用，并自动执行MEX中的代码。codegenclutterSimTracker_kernelarg游戏输入- oclutterSimTracker_kernel;结束

`getTruth`

函数= getTruth(scenario) platpose = platformpose (scenario);%真实信息真理= platPoses (1: end-1);最后一个对象是所有权结束

`clutterSimTracker_kernel`

这个函数定义在一个名为clutterSimTracker_kernel，可在与此脚本相同的工作文件夹中使用。

函数[detectionstracks, ospaMetric] = clutterSimTracker_kernel(Pfas, tgtPoses, insPose, truths, time, systemToRun, systemToReset) assert(numel(Pfas) == 4，'只支持4个错误警报设置。金宝app重写更多持久化变量以添加更多设置'）;持续的Tracker1 tracker2 tracker3 tracker4.．.雷达1雷达2雷达3雷达4.．.．复位1复位2复位3复位4.．.．gospa如果isempty(gospa) || isempty(reset1) || isempty(reset2) || isempty(reset3) || isempty(reset4) gospa = trackgospmetric (“CutoffDistance”, 50);结束如果isempty(tracker1) || isempty(radar1) || isempty(reset1) tracker1 = setupTracker(pfs (1));radar1 = setupRadar (pfa (1));reset1 = 0 (1, 1);结束如果isempty(tracker2) || isempty(radar2) || isempty(reset2) tracker2 = setupTracker(pfs (2));radar2 = setupRadar (pfa (2));reset2 = 0 (1, 1);结束如果isempty(tracker3) || isempty(radar3) || isempty(reset3) tracker3 = setupTracker(pfs (3));radar3 = setupRadar (pfa (3));reset3 = 0 (1, 1);结束如果isempty(tracker4) || isempty(radar4) || isempty(reset4) tracker4 = setupTracker(Pfas(4));radar4 = setupRadar (pfa (4));reset4 = 0 (1, 1);结束开关systemToRun情况下1 detections = radar1(tgtpose, insPose, time);Tracks = tracker1(检测，时间);情况下2 detection = radar2(tgtpose, insPose, time);Tracks = tracker2(检测，时间);情况下3 detection = radar3(tgtpose, insPose, time);Tracks = tracker3(检测，时间);情况下4 detection = radar4(tgtpose, insPose, time);Tracks = tracker4(检测，时间);否则错误(“Idx出界”）;结束[~，gp，~，loc,mT,fT] = gospa(轨道，真理);ospaMetric = [gp loc mT fT];开关systemToReset情况下1 reset1 = 0 (0,1);情况下2 reset2 = 0 (0,1);情况下3 reset3 = 0 (0,1);情况下4 reset4 = 0 (0,1);结束结束函数if (Pfa = Pfa) Pd = 0.9;VCell = 0.0609;细胞体积内联值%;Kc = Pfa / VCell;%杂波密度%的出生率出生率= 0.05*(5 + Pfa*48000)/0.1;变换函数和极限sensorTransformFcn = @ x (x, params);sensorTransformParameters =结构;sensorLimits = bsxfun(@times，[-inf inf]，ones(6,1));将信息组装到trackingSensorConfiguration中配置= trackingSensorConfiguration (“SensorIndex”,1，.．.“IsValidTime”,真的,.．.%更新每一步“DetectionProbability”帕金森病,.．.%可检测状态的检测概率“ClutterDensity”Kc,.．.%杂波密度“SensorTransformFcn”sensorTransformFcn,.．.“SensorTransformParameters”sensorTransformParameters,.．.“SensorLimits”sensorLimits,.．.“MaxNumDetsPerObject”,1，.．.“FilterInitializationFcn”, @initcvgmphd);追踪= trackerPHD (“SensorConfigurations”配置,.．.“出生率”,出生率);结束函数雷达= setupRadar(Pfa) Pd = 0.9;雷达= monostaticRadarSensor (1,.．.“UpdateRate”10.．.“DetectionProbability”帕金森病,.．.“FalseAlarmRate”Pfa,.．.“FieldOfView”, 120年[1],.．.“ScanMode”,“没有扫描”,.．.“DetectionCoordinates”,“场景”,.．.%在场景框架中报告“之内”,真的,.．.%启用INS在场景框架中报告“ReferenceRange”, 300,.．.%短程“ReferenceRCS”10.．.“MaxUnambiguousRange”, 400,.．.%短程“RangeResolution”,1，.．.“AzimuthResolution”,1，.．.“HasFalseAlarms”,真正的);结束