MHT多假设跟踪算法Matlab程序实现:从假设生成到航迹管理
简介这份资源是MHT多假设跟踪算法在Matlab环境下的完整实现程序面向具备Matlab编程基础、正在研究多目标跟踪的开发者与研究生用于解决雷达、视频监控等场景中目标诞生、消失、分裂、合并等不确定情况下的跟踪难题。压缩包共33个文件约45KB以21个.m源码文件为核心辅以4个.c与mexw32、mexglx等混合编程文件另有mat数据、txt说明及html文档覆盖预测、更新、关联、剪枝等模块。已有2037人学习下载。代码围绕卡尔曼滤波预测、概率数据关联、匈牙利与Murty分配、假设生成与剪枝等关键环节展开读者可据此理解MHT从数据预处理、运动建模到关联匹配、分支合并的完整流程并直接复用于自身多目标跟踪项目同时借助可视化脚本分析调试结果。1. MHT算法的Matlab程序从多假设跟踪到可复现的工程实现多假设跟踪Multiple Hypothesis TrackingMHT在工程圈里一直有个尴尬的定位论文里被反复引用真正落地时却常被工程师绕开转而用 JPDA 或 GM-PHD 凑合。原因不复杂——MHT 的假设树会随目标数和杂波密度指数膨胀稍不留神内存就爆了。但它的优势同样硬核在低信噪比、目标交叉、虚警密集的场景下MHT 的航迹起始和维持能力明显优于单假设方法。这篇笔记就围绕「MHT算法的Matlab程序」这个主题把假设生成、剪枝策略、航迹管理到代码落地的完整链路拆开讲。适合已经写过卡尔曼滤波、想进一步啃多目标跟踪的从业者也适合做雷达、红外、视频多目标跟踪方向、需要在 Matlab 里快速验证算法的同学。下面所有代码都是可运行的骨架参数按我实际调过的经验给不追求花哨追求能跑通、能改、能复现。2. MHT算法的核心机制与Matlab实现骨架MHT 的本质不是某一个公式而是一套「延迟决策」的框架。它不急着在每一帧把量测硬分配给航迹而是把所有合理解释保留成假设等后续帧的信息进来后再回溯剪枝。理解这一点后面的代码结构就顺了。2.1 假设树、航迹与量测的三层数据结构MHT 的数据组织通常分三层量测层当前帧的观测点、航迹层已确认或 tentative 的轨迹、假设层一组互斥的航迹-量测关联方案。每一帧新量测到来时对每条已有航迹枚举「关联到某个量测」或「漏检」两种可能再叠加「新目标起始」的分支组合成新的假设集合。在 Matlab 里我一般用结构体数组加嵌套 cell 来表达而不是硬写类。原因是调试时disp一眼能看清层级改字段也方便。核心字段包括tracks每条航迹的状态均值、协方差、航迹 ID、生命周期计数、hypotheses每个假设包含的关联矩阵和累计似然、meas当前帧量测矩阵每行一个[x y]或[x y z]。% 初始化MHT全局结构 mht.tracks struct(id, {}, state, {}, P, {}, life, {}, score, {}); mht.hypotheses struct(assoc, {}, logLik, {}, parent, {}); mht.meas []; % 当前帧量测 Nx2 mht.nextId 1; % 航迹ID分配器 mht.gateThresh 9.21; % 卡方门限2自由度95%对应9.21 mht.pruneThresh -20; % 假设对数似然剪枝门限 mht.confirmHits 3; % 连续命中3次确认航迹 mht.deleteMiss 4; % 连续漏检4次删除航迹这段结构体定义是整个程序的「黑匣子」入口。gateThresh用卡方分布门限做波门2 自由度 95% 置信对应 9.213 自由度对应 11.34按你的量测维度选。pruneThresh是对数似然剪枝线设得太松假设爆炸太紧会误删正确关联我一般从 -20 起步观察假设数量再调。confirmHits和deleteMiss是航迹生命周期管理M 取 3、N 取 4 是工程上比较稳的组合虚警多时把 confirmHits 提到 4。2.2 假设生成枚举关联与似然计算假设生成是 MHT 最耗算力的环节。对每条航迹先用马氏距离做波门筛选只保留落在门内的量测再枚举关联。似然计算用高斯量测模型漏检概率和虚警密度作为先验参与打分。function [newHyps] generateHypotheses(mht, Pd, lambda) % Pd: 检测概率, lambda: 虚警密度(每帧每单位面积) newHyps struct(assoc, {}, logLik, {}, parent, {}); for h 1:numel(mht.hypotheses) baseAssoc mht.hypotheses(h).assoc; baseLik mht.hypotheses(h).logLik; % 对每条航迹独立枚举漏检 或 关联到门内某个量测 branchList cell(1, numel(mht.tracks)); for t 1:numel(mht.tracks) trk mht.tracks(t); opts {{miss, 0}}; % 漏检分支 for m 1:size(mht.meas,1) d mht.meas(m,:) - trk.state(1:2); S trk.P(1:2,1:2) eye(2)*0.5; % 量测噪声R md d / S * d; if md mht.gateThresh opts{end1} {meas, m, md}; %#okAGROW end end branchList{t} opts; end % 笛卡尔积组合所有航迹的分支 combos cartesianBranches(branchList); for c 1:numel(combos) combo combos{c}; logL baseLik; assoc baseAssoc; usedMeas []; valid true; for t 1:numel(combo) if strcmp(combo{t}{1}, miss) logL logL log(1 - Pd); else m combo{t}{2}; if ismember(m, usedMeas) % 一个量测不能分给两条航迹 valid false; break; end usedMeas(end1) m; %#okAGROW md combo{t}{3}; logL logL log(Pd) - 0.5*md - log(2*pi*det(S)^0.5); assoc(t) m; end end if ~valid, continue; end % 未关联量测视为新目标起始或虚警 nNew size(mht.meas,1) - numel(usedMeas); logL logL nNew*log(lambda); newHyps(end1) struct(assoc, assoc, logLik, logL, parent, h); %#okAGROW end end end逻辑上分三步先对每条航迹建候选分支漏检 门内量测再做笛卡尔积组合最后逐组合校验「一个量测只能分配给一条航迹」并累加对数似然。Pd是检测概率雷达场景常取 0.9视频检测器取 0.7~0.85。lambda是虚警密度按你场景里每帧平均虚警数除以视场面积估。这里用对数似然而非概率是为了避免连乘下溢剪枝时直接比大小即可。注意S里我加了量测噪声R实际项目里 R 要从传感器手册或标定数据来别拍脑袋。2.3 剪枝与航迹管理让假设树不爆炸假设数量是 MHT 的命门。三条航迹、每条门内两个量测组合数就是 3 的 3 次方 27五条航迹就上百。必须剪枝。function mht pruneAndManage(mht) % 1) 按对数似然排序保留前K个假设 [~, idx] sort([mht.hypotheses.logLik], descend); K 50; % 假设上限按内存调 idx idx(1:min(K, numel(idx))); mht.hypotheses mht.hypotheses(idx); % 2) 删除低于剪枝门限的假设 keep [mht.hypotheses.logLik] mht.pruneThresh; mht.hypotheses mht.hypotheses(keep); if isempty(mht.hypotheses) mht.hypotheses struct(assoc, [], logLik, 0, parent, 0); end % 3) 航迹生命周期更新以最优假设为准 best mht.hypotheses(1); for t 1:numel(mht.tracks) if best.assoc(t) 0 mht.tracks(t).life(1) mht.tracks(t).life(1) 1; % 命中 mht.tracks(t).life(2) 0; else mht.tracks(t).life(2) mht.tracks(t).life(2) 1; % 漏检 end % 确认与删除 if mht.tracks(t).life(1) mht.confirmHits mht.tracks(t).score 1; % 确认 end if mht.tracks(t).life(2) mht.deleteMiss mht.tracks(t).score -1; % 待删除 end end mht.tracks mht.tracks([mht.tracks.score] 0); end剪枝分两层先按似然保留前 K 个假设K 取 50 是内存和精度的折中目标多时提到 100再砍掉低于门限的。航迹管理用「连续命中/连续漏检」计数比单纯累计分数更抗虚警。这里有个血泪经验life字段我用两元素向量存命中数和漏检数别用单个计数器否则漏检后命中无法正确清零航迹会「赖着不死」。删除航迹时同步清理假设里的关联索引否则后续帧会访问越界。3. 用Matlab跑通MHT从仿真数据到航迹输出光有骨架不够得让它跑起来看到航迹。这一章用一个匀速直线运动的仿真场景把数据生成、主循环、卡尔曼更新串成完整可执行流程。3.1 仿真场景搭建匀速目标加杂波我一般先写一个数据生成脚本目标数、帧数、杂波密度都可调方便后面做消融。function [measSeq, truth] genScenario(nFrame, nTarget, lambda, Pd) % 生成匀速直线目标 均匀杂波 F [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1]; % 状态转移 H [1 0 0 0; 0 1 0 0]; % 量测矩阵 Q diag([0.01 0.01 0.1 0.1]); % 过程噪声 R diag([0.5 0.5]); % 量测噪声 truth cell(1, nTarget); for i 1:nTarget x [rand*50; rand*50; (rand-0.5)*4; (rand-0.5)*4]; truth{i} zeros(4, nFrame); for k 1:nFrame truth{i}(:,k) x; x F*x mvnrnd(zeros(4,1), Q); end end measSeq cell(1, nFrame); for k 1:nFrame pts []; for i 1:nTarget if rand Pd z H*truth{i}(:,k) mvnrnd(zeros(2,1), R); pts [pts; z]; %#okAGROW end end nClutter poissrnd(lambda*100); % 视场100x100 clutter [rand(nClutter,1)*100, rand(nClutter,1)*100]; measSeq{k} [pts; clutter]; end endF和H是标准匀速模型Q控制目标机动程度R是量测噪声。杂波用泊松分布生成lambda*100是每帧期望虚警数视场 100×100。这里Pd直接决定目标点是否出现方便你观察检测概率对 MHT 的影响。注意mvnrnd每次调用有开销帧数上千时建议预生成随机数否则仿真本身比算法还慢。3.2 主循环预测、关联、更新一条龙主循环把前面几块拼起来顺序是「预测 → 生成假设 → 剪枝 → 用最优假设做卡尔曼更新」。function runMHT() [measSeq, truth] genScenario(60, 3, 0.9, 0.5); mht initMHT(); F [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1]; H [1 0 0 0; 0 1 0 0]; Q diag([0.01 0.01 0.1 0.1]); R diag([0.5 0.5]); for k 1:numel(measSeq) mht.meas measSeq{k}; % 1) 预测 for t 1:numel(mht.tracks) mht.tracks(t).state F * mht.tracks(t).state; mht.tracks(t).P F * mht.tracks(t).P * F Q; end % 2) 假设生成 剪枝 mht.hypotheses generateHypotheses(mht, 0.9, 0.5); mht pruneAndManage(mht); % 3) 用最优假设做卡尔曼更新 best mht.hypotheses(1); for t 1:numel(mht.tracks) if best.assoc(t) 0 z mht.meas(best.assoc(t), :); S H * mht.tracks(t).P * H R; K mht.tracks(t).P * H / S; mht.tracks(t).state mht.tracks(t).state K*(z - H*mht.tracks(t).state); mht.tracks(t).P (eye(4) - K*H) * mht.tracks(t).P; end end % 4) 新目标起始未关联量测建 tentative 航迹 used best.assoc(best.assoc 0); for m 1:size(mht.meas,1) if ~ismember(m, used) nt.state [mht.meas(m,:); 0; 0]; nt.P diag([1 1 10 10]); nt.id mht.nextId; mht.nextId mht.nextId 1; nt.life [1 0]; nt.score 0; mht.tracks(end1) nt; %#okAGROW end end % 5) 记录确认航迹用于画图 confirmed mht.tracks([mht.tracks.score] 1); fprintf(Frame %d: %d confirmed tracks\n, k, numel(confirmed)); end end顺序不能乱预测必须在假设生成前否则波门用的是旧状态。卡尔曼更新只对最优假设执行这是 MHT「延迟决策」的体现——其他假设保留但不更新状态等后续帧翻盘。新目标起始用未关联量测初始速度设 0、速度方差设大这里 10让滤波器自己收敛。fprintf那行是调试用的看确认航迹数是否稳定如果一直涨说明虚警在起航迹得提高confirmHits或降lambda。3.3 关键参数怎么调一张对照表参数调不好MHT 要么丢目标要么炸内存。下表是我在不同场景下反复试出来的经验区间。参数含义稀疏杂波密集杂波调大后果gateThresh波门卡方门限9.2111.34门内量测多假设膨胀pruneThresh对数似然剪枝线-15-25保留假设多精度升但变慢K假设数上限30100内存和耗时线性涨confirmHits确认所需命中34起航迹慢抗虚警强deleteMiss删除所需漏检43航迹存活久易残留调参顺序建议先定gateThresh保证正确量测落在门内用真值算马氏距离分布再调pruneThresh和K控制假设规模最后用confirmHits/deleteMiss平衡起航速度和虚警。别一上来就调剪枝门限波门不对后面全白搭。4. MHT程序落地时的避坑与排查清单MHT 代码跑不通八成不是算法理解错而是几个工程细节翻车。下面五条是我踩过的坑按「现象 → 原因 → 解决」写。现象假设数量每帧翻倍几帧后内存爆掉。原因笛卡尔积组合时没做「一个量测只能分配给一条航迹」的互斥校验或者剪枝在假设生成之后才做中间态已经爆炸。解决在组合循环里用usedMeas实时校验一旦冲突立即break同时把剪枝上限 K 设小一点先跑通再逐步放大。现象航迹在目标真实位置附近反复起起灭灭ID 一直变。原因confirmHits设得太低虚警也能凑够命中数或者新目标起始时初始协方差太小滤波器把杂波当目标锁死。解决confirmHits提到 4新航迹初始速度方差放大到 10~20让它在几帧内被真实量测拉过去拉不过去自然消亡。现象目标交叉时航迹互换两条轨迹「粘」在一起。原因MHT 的最优假设选择只看累计似然交叉瞬间两个关联方案似然接近选错后没有回溯机制。解决保留多个假设更久提高 K 和放宽 pruneThresh并在航迹更新时对交叉区域的目标用更大的过程噪声 Q让滤波器对机动更宽容。这是 MHT 相对 JPDA 的优势所在别浪费。现象程序越跑越慢后期每帧耗时是前期的十倍。原因mht.tracks和mht.hypotheses用end1动态增长Matlab 反复重分配内存或者删除航迹后没清理假设里的索引导致无效计算。解决预分配结构体数组删除时用逻辑索引一次性重建假设里的assoc索引在航迹删除后要同步重映射否则会访问到已删除航迹。现象仿真里航迹完美换真实数据就崩。原因仿真量测噪声是各向同性的高斯真实传感器噪声有偏、有相关性R 矩阵设错导致波门失效。解决用真实数据先做量测残差统计估出实际 R 再填进去如果噪声非高斯波门门限要放宽或者改用椭圆波门加自适应门限。提示调试 MHT 时先把Pd设成 1、lambda设成 0确认纯关联逻辑正确再逐步加杂波和漏检。这样出问题能快速定位是关联错还是鲁棒性不够。5. 让MHT跑得更稳假设回溯与航迹融合的进阶技巧基础版跑通后真正拉开差距的是两个进阶点假设回溯和航迹融合。前者解决「当前帧选错关联、后续帧能否纠正」后者解决「多条假设对应同一目标时如何输出唯一航迹」。假设回溯的做法是给每个假设记录parent指针前面代码里已经留了字段当某个假设的累计似然在后续帧反超时沿着 parent 链回溯把中间帧的关联决策一并修正。实现上不需要真的改历史状态只需在输出航迹时按最优假设的祖先链重新走一遍卡尔曼更新。代价是每个假设要存关联历史内存换精度。我一般只在目标数少于 10 的场景开回溯目标多了内存吃不消。航迹融合更实用。MHT 的多个假设里同一目标可能对应多条 tentative 航迹直接按最优假设输出会漏掉一些真实目标。常见做法是对确认航迹做一次「相似度合并」状态马氏距离小于门限、且速度方向夹角小于 30 度的两条航迹按似然加权融合成一条。function tracks mergeTracks(tracks, gate) % 合并状态接近的确认航迹gate为马氏距离门限 i 1; while i numel(tracks) j i 1; while j numel(tracks) d tracks(i).state(1:2) - tracks(j).state(1:2); S tracks(i).P(1:2,1:2) tracks(j).P(1:2,1:2); if d / S * d gate w tracks(i).life(1) / (tracks(i).life(1) tracks(j).life(1)); tracks(i).state w*tracks(i).state (1-w)*tracks(j).state; tracks(i).P w^2*tracks(i).P (1-w)^2*tracks(j).P; tracks(j) []; else j j 1; end end i i 1; end end融合权重用命中次数占比命中多的航迹话语权大。gate取 4~6 比较稳太小合并不了太大把邻近的真实双目标误合。这个函数放在剪枝之后、输出之前调用能明显减少「一个目标两条航迹」的冗余。验证 MHT 效果别只看航迹图要算两个量化指标一是 OSPA 距离最优子模式分配衡量估计航迹与真值的整体偏差二是航迹 ID 切换次数反映关联稳定性。OSPA 在 Matlab 里可以自己实现核心是枚举估计与真值的分配组合取最小代价目标少时直接暴力枚举即可。我习惯每帧记录确认航迹状态和真值跑完算平均 OSPA 和总切换次数调参时盯着这两个数比肉眼看图靠谱得多。最后说个习惯MHT 的参数没有普适最优换场景就得重调。我一般把gateThresh、pruneThresh、K三个参数做成脚本顶部的变量跑之前先做一轮网格搜索用 OSPA 选最优组合再固化到配置里。别指望一套参数打天下那是玄学不是工程。希望帮到你。本文还有配套的精品资源点击获取