简介基于随机集理论的多目标跟踪算法Matlab工具包涵盖PHD、CPHD、CBMeMBer等经典滤波器实现适合研究雷达、视频监控等领域多目标跟踪的算法工程师与研究生。压缩包共206个文件其中205个为.m脚本包含滤波主程序、数据关联如Hungarian算法与可视化模块另有1个说明文档整体仅166KB结构精简便于快速移植。已有767人学习下载适合作为从理论到代码的桥梁。通过run_filter等示例脚本可掌握随机集滤波器的更新、预测与航迹提取流程快速复现文献结果并开展对比实验工具包源于VoBN的工作代码风格清晰注释与模块划分较规范结合标签与描述可看出其典型应用场景。整体上这是一份实践价值较高的多目标跟踪算法代码集能帮助读者节省编码时间聚焦算法核心逻辑。1. 随机集理论为什么能统一多目标跟踪从 PHD 到 CBMeMBer雷达、无人机集群和视频监控里目标数量随时会变检测器还会漏报、误报传统多目标跟踪必须先做数据关联目标一多、杂波一密分配组合数就会爆炸。随机集理论把每个目标状态当成集合元素把量测也当成集合用贝叶斯递推直接估计集合后验绕开了显式关联这一环。PHD 滤波器只传递强度函数的一阶矩CPHD 额外估计目标基数分布CBMeMBer 则用多伯努利参数近似并做势平衡。这套 matlab 工具包把三种算法放在同一套仿真框架下入口是 run_filter.m关联部分用 Hungarian.m非常适合做随机集滤波算法对比和科研复现。对需要快速上手的工程师来说跑通 demo 再改参数比从论文公式推代码要高效得多。2. 随机集建模与工具包文件结构PHD/CPHD/CBMeMBer 的三条技术路线2.1 从 RFS 到概率假设密度三种滤波器在近似方式上的区别多目标跟踪的难点在于状态空间和量测空间的维数都在随时间变化。随机有限集把 k 时刻所有单目标状态写成集合 (X_k{x_{k,1},...,x_{k,M_k}})量测写成集合 (Z_k{z_{k,1},...,z_{k,N_k}})目标数量变化就转化为集合基数的变化。基于这个建模贝叶斯递推可以在集合积分意义下展开但完整后验分布仍然无法直接计算于是需要做近似。PHD 滤波器用泊松随机有限集来近似后验只传播强度函数 v(x)也就是后验的一阶矩。对 v(x) 在状态空间上积分可以得到目标数量的期望。泊松假设意味着目标数分布的方差等于均值杂波密度较高时这个假设会把方差估计偏大导致滤波器目标数在部分帧里抖动明显。我第一次用 PHD 跑密集杂波场景时真实目标数是 5滤波器估计值却在 3 到 9 之间来回跳后来对照 CPHD 才确认是泊松近似的固有偏差。CPHD 的改进是额外维护一个基数分布 p_k(n)表示 k 时刻目标数为 n 的概率分布。更新时联合估计基数分布和强度函数目标数方差不再被强制等于均值所以杂波不均匀的场景下精度更高但每帧都要做基数分布的卷积运算计算量比 PHD 高一个量级。CBMeMBer 则改用多伯努利随机有限集近似后验每个候选航迹都有存在概率 r 和空间分布 p(x)更新后通过势平衡修正目标数偏差。它的优势是不依赖泊松假设目标数大幅跳变时表现更稳实际代码里需要维护伯努利分量数组剪枝和合并策略对最终效果影响很大。为了快速对比选型我习惯用下面这张表梳理核心差异滤波器后验近似形式关键假设计算复杂度常见短板PHD泊松 RFS仅强度函数目标数服从泊松分布低O(NM)目标数方差被高估CPHD强度函数 基数分布基数分布独立于状态中高涉及基数卷积检测概率不均时复杂度高CBMeMBer多伯努利 RFS各航迹存在概率独立中高分量数随目标增长需要精细剪枝合并参数这张表只覆盖理论假设。实际工程中如果量测是非线性高斯模型三种滤波器都要配合粒子滤波或无迹变换工具包里的 matlab 代码默认走线性高斯路径换场景时需要同步改模型矩阵。2.2 工具包文件构成run_filter.m 作为统一入口解压 zip 后看到大量重复的 run_filter.m这是 VoBN 系工具包的典型组织方式不同滤波器放在不同子目录里每个目录下各有一个 run_filter.m 作为该算法的仿真入口Hungarian.m 则放在公共目录里供多个滤波器调用。由于 MATLAB 按 path 顺序搜索文件直接把所有子目录一次性 addpath运行时可能调用到错误版本这点在后续跑 demo 时需要特别留意。run_filter.m 的逻辑一般分为四段定义运动模型和量测模型初始化滤波器参数结构体进入时间循环递推最后绘图并输出误差指标。下面这段代码是这类脚本最常见的骨架我以 PHD 为例做了简化% run_filter.m 主入口骨架以 PHD 为例 simulation_length 100; % 仿真帧数 state_dim 4; % [x vx y vy] birth_model struct(mean, [0;0;0;0], ... cov, eye(state_dim)*100, ... weight, 1.0); % 新生目标强度 % 滤波器参数结构体 tracker struct(lambda_c, 10, ... % 杂波泊松均值 P_D, 0.98, ... % 检测概率 P_G, 0.999, ... % 波门概率 model, model, ... % 匀速运动模型 birth, birth_model); for k 1:simulation_length % 生成真实多目标状态 X_true 和量测 Z [X_true, Z] simulate_scene(tracker, k); % 调用 PHD 滤波器主体 [tracker, est_state] phd_filter(tracker, Z, k); % 绘制每一帧结果并累积 OSPA 误差 plot_frame(X_true, est_state, k); end代码里的lambda_c是每帧量测中杂波点数的期望值P_D是传感器对单目标的检测概率P_G决定波门覆盖量测的百分比。simulate_scene内部包含目标出生、死亡、匀速运动和量测噪声生成这些在工具包里对应单独的 m 文件。如果直接复制 run_filter.m 到新场景必须同步检查model结构体里的状态转移矩阵和过程噪声协方差否则可能在头几帧就发散。之前我修改 P_D 时忘了换运动模型CPHD 的目标数估计在前 20 帧内跳到了真实值的四倍排查半天才发现是模型矩阵和场景不匹配。3. 在 MATLAB 中跑通多目标跟踪 demo运行环境与参数调整3.1 运行环境与第一个 demo 的启动方式这套工具包对 MATLAB 版本要求不算高R2017b 以上通常就能直接运行不需要额外的 Simulink 或深度学习工具箱。代码大量使用矩阵操作和 cell 数组只要基础环境完整即可。唯一要注意的是解压路径不能带中文或空格否则脚本里的相对路径加载会失败报错也经常是含糊的“文件不存在”。启动流程分四步第一步把压缩包解压到纯英文目录第二步在 MATLAB 中 cd 到目标子目录第三步执行clear; close all; clc;清空工作区第四步直接运行run_filter。正常时会弹出逐帧更新的轨迹图同时命令行打印当前帧的目标数和 OSPA 误差值。为了复现相同实验我通常先固定随机种子% 固定随机数种子便于复现实验 rng(2024, twister); run_filter固定种子是因为工具包中的目标出生位置、量测噪声和杂波位置都是随机生成的。不固定种子同一组参数跑两次 OSPA 曲线会有明显差异调参时很难判断是参数生效还是随机波动。twister是 Mersenne Twister 随机数生成器从 R2011a 开始就是 MATLAB 默认算法。后续做 Monte Carlo 实验时外层循环里每轮换一次种子最后对误差取平均才能得到平滑的性能曲线。3.2 关键参数配置杂波、检测概率与新生目标调整实验场景不需要大改滤波器主体多数情况下只需要修改 run_filter.m 开头的结构体字段。这里是我整理过的几组关键参数及调参观察点参数典型值作用调参时观察的指标lambda_c515每帧杂波数目的泊松均值虚警率、航迹碎片数P_D0.90.99目标被正确检测的概率航迹连续性、目标数误差P_G0.999量测波门概率控制关联候选集OSPA 距离突变时刻birth_model.covdiag([100,10,100,10])新生目标状态协方差航迹起批速度和初次误跟以lambda_c为例把它从 10 改成 50每帧杂波点增加到原来的五倍。PHD 滤波器会把部分杂波当成低权重新生目标目标数估计波动明显变大CPHD 因为基数分布修正目标数均值能稳定很多但基数分布的收敛帧数会拉长CBMeMBer 则必须同步提高剪枝阈值否则伯努利分量数量膨胀运行耗时直线上升。调整P_D的影响更微妙降到 0.8 时漏检概率提升PHD 强度函数在漏检帧被整体衰减目标数积分会瞬间掉一截。缓解手段是适当调高存活概率P_S或者给强度函数设置一个最小下限避免航迹在短暂漏检后彻底丢失。新生目标协方差birth_model.cov不能取得太小。若位置方差设置小于量测噪声方差新目标的强度会被滤波器当成杂波过滤掉航迹起批会延迟好几帧。我一般取场景中位置噪声方差的 2 到 3 倍速度项根据目标最大机动速度估算。这样在目标进入雷达范围的初始帧滤波器就能迅速建立新航迹不会等目标完全进入波门后才反应。4. Hungarian.m 在数据关联中的角色与实现4.1 Hungarian.m 的函数接口与调用方式随机集滤波器在更新阶段避开了显式关联但输出连续航迹时仍然要把不同帧的估计状态对应起来。Hungarian.m 在这个工具包里就是用来做线性指派的最小代价匹配。它接收一个 m 行 n 列的代价矩阵返回最优关联矩阵和最小总代价算法复杂度为 O(n^3)n 取 m 和 n 中的较大值。在多目标跟踪里行通常对应当前帧提取的估计状态列对应上一帧保留下来的航迹预测状态。代价矩阵里的元素一般是马氏距离因为欧氏距离没有考虑状态协方差。下面是一段调用 Hungarian.m 的典型代码% 构造代价矩阵行为航迹列为当前帧量测 costMat zeros(numTracks, numMeasurements); for i 1:numTracks for j 1:numMeasurements diff meas(:,j) - track_pred(:,i); costMat(i,j) diff / track_cov * diff; end end % 波门限制超过门限的候选代价置为无穷 costMat(costMat GATE_THRESHOLD) Inf; % 调用 Hungarian 算法返回最优关联矩阵 [assignment, totalCost] Hungarian(costMat);其中track_pred是上一帧航迹状态推进到当前时刻的预测值track_cov是预测协方差diff是量测与预测的差值。GATE_THRESHOLD通常取chi2inv(0.99, state_dim)也就是卡方分布的 0.99 分位数用来控制允许进入匹配候选集的量测比例。assignment是 numTracks 行 numMeasurements 列的 0-1 矩阵约束是每行每列至多一个 1totalCost则是最优匹配下的马氏距离总和可作为关联质量的度量。有一点需要留意匈牙利算法要求代价矩阵非负工具包内部会把矩形矩阵补成方阵补上的虚行虚列代价设为大数。如果杂波点远多于真实航迹数方阵规模会由量测数决定计算耗时上涨很快。之前处理过量测 200、航迹 10 的场景单次 Hungarian 耗时从 0.2 毫秒涨到 3 毫秒100 帧仿真就多出近 0.3 秒。对实时系统来说应加门控预筛而不是直接丢给匈牙利算法暴力求解。4.2 关联代价的构造与门控策略匈牙利算法不关心代价来源但代价矩阵定义直接影响关联质量。若直接用位置欧氏距离状态里位置和速度量纲不同速度分量几乎不参与代价计算关联很容易把位置接近但速度差异大的目标错配。马氏距离通过协方差矩阵归一化能同时考虑位置和速度不确定性。具体做法是计算量测预测协方差 S然后用diff * inv(S) * diff作为二次型距离。门控策略是在构造代价矩阵前筛掉距离上明显不可能的配对。椭圆波门的判定条件就是马氏距离小于卡方分布阈值。把这段逻辑向量化比双重 for 循环高效% 向量化计算马氏距离避免双重循环 S_inv inv(track_cov); % d2 为 numMeasurements x numTracks 的二次型距离矩阵 d2 zeros(nMeas, nTracks); for j 1:nTracks delta meas - track_pred(:,j); % 每一列是一个量测 d2(:,j) sum(delta .* (S_inv * delta), 1); end % 波门将大于阈值的项标记为不可行 d2(d2 chi2inv(0.99, state_dim)) Inf; [assignment, cost] Hungarian(d2);这里的track_cov是当前航迹的量测预测协方差state_dim由量测函数决定。取 0.99 分位数意味着在正确量测服从高斯分布的前提下最多允许百分之一的正确量测被门限排除。门限取太大杂波大量进入候选集匈牙利算法要处理很多无效列门限取太小真实量测被排除后续航迹更新会失去观测支持。常见的错误是代价矩阵全为 Inf 时部分 Hungarian.m 实现会返回空关联矩阵后续更新代码如果没有判空就直接索引会崩溃。我的做法是在调用前检查all(isinf(costMat(:)))是则跳过该帧关联保留航迹预测状态。下面给出 Hungarian.m 接口边界条件的速查表项目说明costMat尺寸m 行 n 列m 为航迹数n 为量测数内部补边虚行虚列代价设为大数assignment约束每行每列至多一个 1复杂度O(n^3)n 为 max(m,n)常见错误全 Inf 矩阵导致空分配5. 验证多目标跟踪效果的三个进阶技巧OSPA 距离与滤波一致性5.1 用 OSPA 距离对比 PHD/CPHD/CBMeMBer多目标误差不能简单对每个目标单独求误差再求和因为滤波器估计的目标数和真实目标数经常不一致。OSPA 距离是当前最常用的集合距离关键参数是截断距离 c 和阶次 p。c 控制单个目标误差的最大惩罚上限p 控制误差范数的敏感度通常取 c100、p2。若工具包没有内置 OSPA 计算函数可以按下面框架自己写function d ospa_metric(X, Y, c, p) % X, Y 分别为真实和估计状态集合每行是一个目标状态 m size(X,1); n size(Y,1); if m n [X,Y] deal(Y,X); [m,n] deal(n,m); end D zeros(m,n); for i 1:m for j 1:n D(i,j) min(norm(X(i,1:2)-Y(j,1:2)), c); end end [assign, ~] Hungarian(D); cost 0; for k 1:m idx assign(k,:) 1; if any(idx) cost cost D(k, idx)^p; else cost cost c^p; end end d (cost / n)^(1/p); end这个函数里再次用到了 Hungarian目的是寻找真实目标和估计目标间的最优距离配对。如果 m n说明估计目标不足未匹配的 n-m 个目标每个直接贡献 c^p。最终结果除以 n 后开 p 次方这样目标数误差和位置误差被折算进同一个标量。用这个指标对比三种滤波器时要固定同一组量测数据否则目标出生和杂波位置不一致差异会淹没算法本身的差距。5.2 参数敏感性扫描的通用做法单次仿真只能说明一个场景下的表现更有用的是扫描关键参数观察算法在什么条件下开始失效。常见做法是在 run_filter 外层循环lambda_c和P_D每组参数固定随机种子跑一次记录 OSPA 均值% 外层扫描脚本自动保存 OSPA 结果 lambda_range [5 10 20 50]; pd_range [0.8 0.9 0.98]; ospa_table zeros(length(lambda_range), length(pd_range)); for i 1:length(lambda_range) for j 1:length(pd_range) tracker.lambda_c lambda_range(i); tracker.P_D pd_range(j); rng(100 i*10 j); % 每组不同种子 run_filter; % 内部自动计算 ospa ospa_table(i,j) mean(ospa_over_time); end end每组参数的随机种子必须不同否则不同配置下的量测噪声和杂波位置完全相同只是改了滤波参数测出来的差异实际上是同一组噪声下的人为结果会高估参数灵敏度。扫描完成后用imagesc(ospa_table)画热图可以快速看出算法在哪个区域性能突然恶化。除了 OSPA还可以把 Hungarian.m 的代价矩阵输出保存到文件复盘关联错误发生的具体时帧。这个调试矩阵不参与主流程却能让整个工具包从黑盒变成可观测的仿真环境排查滤波发散问题时比只看误差曲线直接得多。本文还有配套的精品资源点击获取