上次帮一个师弟调雷达目标跟踪的仿真,他卡在了一个很典型的问题上:状态方程明明是线性的,量测方程里只要出现一个atan2,他就不敢套卡尔曼滤波了。这个问题问得非常有代表性。在滤波跟踪这个领域,很多人能把EKF、UKF、PF三个缩写背得滚瓜烂熟,也知道一个线性化、一个无迹变换、一个蒙特卡洛,但等到真要把它们放到同一个量测非线性模型下跑一遍、用同一把尺子量出高下的时候,往往不知道该从哪儿下手。这篇博文就想把这层窗户纸捅破:我用一个极坐标量测的匀速运动目标模型做仿真场景,把扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)和粒子滤波(PF)放到同一个赛场上,从滤波精度、一致性和计算耗时三个维度做了对比,并把Matlab实现的关键代码和调参经验一并整理出来。无论你是刚接触非线性滤波的研究生,还是要在工程里做目标跟踪选型的开发人员,这篇文章都能帮你少走不少弯路。
1. 为什么题目落在“量测非线性”这个要害上
先说说为什么这个题目把重点放在“量测非线性”上,而不是状态方程非线性。很多教材讲EKF和UKF的区别时,喜欢把状态方程和量测方程混在一起说,导致读者分不清到底哪一步在起作用。实际上,在目标跟踪最常见的场景里,目标运动模型通常建立在地球直角坐标系下,匀速、匀加速转弯这些模型本质上都是线性的,或者可以写成线性的状态转移矩阵。但传感器报告的数据往往是极坐标或者球坐标下的量测——雷达给的是距离、方位角、俯仰角,你用来做跟踪的状态却是X、Y、Z和对应的速度分量。从量测到状态,这中间天然隔着一个非线性映射。
1.1 线性卡尔曼的舒适区与边界
经典卡尔曼滤波有两个硬性前提:系统方程和量测方程都是线性的,噪声服从高斯分布。在这两个前提下,卡尔曼滤波给出的是最小均方误差意义下的最优估计,这个“最优”是有严格数学证明的。可是一旦量测方程里混进了sin、cos、atan2、平方根这些非线性函数,高斯分布经过非线性变换之后就不再是高斯分布了,卡尔曼滤波的正交投影理论体系就崩塌了。你硬用线性框架去套非线性问题,结果就是协方差矩阵根本表达不了真实分布的形状,滤波精度急剧下降,严重的时候直接发散。
1.2 非线性量测在现实中的三个典型来源
做工程的人一定要建立这个直觉:量测非线性不是教科书里的抽象概念,它无处不在。最常见的来源有三个。第一是坐标转换,就是上面说的雷达极坐标量测转直角坐标状态,这是最典型也最容易理解的。第二是角度测量,比如测向交叉定位,量测模型里会出现方位角的正切关系,这种非线性函数在某些角度区间上雅可比矩阵变化非常剧烈,属于比较难处理的场景。第三是带有物理约束的量测模型,比如多普勒雷达的量测方程里状态和量测之间是速度投影关系,再比如某些传感器输出是距离的平方或者信号强度倒数的对数,这些模型天然非线性的程度都很高。
1.3 三种算法分别代表哪三条技术路线
EKF、UKF、PF面对同一个非线性量测问题,走的路线完全不同。EKF的思路是“局部线性化”,用一阶泰勒展开把非线性函数在估计点附近掰直,掰直之后继续沿用卡尔曼滤波的框架。UKF的思路是“确定性采样”,既然一个高斯分布经过非线性变换后形状会变,那就用一组精心构造的Sigma点去捕捉这个分布的关键统计量,通过非线性函数之后再加权合成新的均值和协方差。PF的思路是“蒙特卡洛逼近”,干脆不用高斯分布去近似真实分布,而是用一簇带有权重的随机粒子去直接拟合任意形状的概率密度。这三条路线分别对应着对非线性问题的三种不同妥协方式,也决定了它们在精度、计算量和适用范围上的根本差异。理解了这一层,后面看仿真结果的时候就不会只停留在“谁好谁差”的表面结论上。
2. 三种滤波算法的底层差异:泰勒展开、Sigma点与粒子群
在动手写代码之前,有必要把三种算法的核心机制掰开揉碎讲清楚。很多人跑仿真的时候只是照着开源代码把结果跑出来,却说不清每个步骤在做什么,这样一旦结果不符合预期,根本不知道去哪里排查。我不打算堆一堆公式推导,而是用尽量直白的方式说清楚每一行代码背后的物理意义。
2.1 EKF:把非线性函数在均值点做一阶线性化
EKF的核心假设是:虽然量测函数h(x)是非线性的,但如果我们已经对状态有一个估计,真实状态大概率在这个估计附近的一个小邻域内,那么在这个邻域里用切线去逼近曲线,误差是可以接受的。
具体操作就是对量测函数求雅可比矩阵:
H = [px/r, py/r, 0, 0; -py/r^2, px/r^2, 0, 0];这个矩阵的含义是:当状态变量的每个分量发生微小变化时,量测值(距离和角度)分别会变化多少。它本质上是在用线性关系去近似非线性关系,把h(x)在高斯分布均值点处的切平面当作真实的量测曲面。
EKF的优点是计算量小,工程实现简单,在非线性程度不高的场景下精度尚可。缺点是它只取了一阶项,把高阶信息全丢掉了。如果量测非线性很强,比如目标距离很近时角度变化非常剧烈,或者目标的方位角接近正负90度时atan2的导数变化很快,EKF的线性化误差会被显著放大,协方差的传递就会失真,最终导致滤波发散。此外,EKF要手动推导雅可比矩阵,模型稍微改一改就要重新算导数,工程维护成本也不低。
2.2 UKF:用一组Sigma点代替整个概率分布
UKF的思想比EKF优雅得多。它不去算雅可比,而是回答一个问题:如果x满足均值为x_hat、协方差为P的高斯分布,经过非线性函数h(x)之后,得到的z的均值和协方差大概是多少?
这个问题的难点在于,高斯分布经过非线性变换后不再是高斯分布,解析解很难求。UKF的巧妙之处在于:我不用管整个分布的形状,只需要选一组特殊的点,让它们的均值和协方差恰好等于x的均值和协方差,然后把这组点分别代入非线性函数,得到一组变换后的点,再对这组变换后的点计算加权均值和加权协方差。这一套操作称为无迹变换(Unscented Transform)。
Sigma点的生成方式如下:
lambda = alpha^2 * (n + kappa) - n; P_sqrt = chol((n + lambda) * P, 'lower'); X(:,1) = x; for i = 1:n X(:,i+1) = x + P_sqrt(:,i); X(:,i+n+1) = x - P_sqrt(:,i); end这里alpha控制Sigma点距离均值的远近,通常取一个比较小的正数比如1e-3;kappa是次级缩放因子,高斯分布下取0;beta用于引入先验分布信息,高斯分布下取2最优。这三个参数直接决定了UT变换的精度,也是后面要重点说的调参对象。
UKF的精度至少能捕捉到非线性函数二阶泰勒展开的信息,部分场景下能达到三阶精度,并且完全不需要推导雅可比矩阵。它的计算量大约是EKF的2到3倍,对于现代计算机来说几乎可以忽略不计。在大多数量测非线性场景下,UKF用EKF相同的计算成本换来了更好的精度和鲁棒性,这也是它在工程里越来越受青睐的原因。
2.3 PF:用带权重的粒子群做直接逼近
粒子滤波的思想跟前面两者完全不一样。EKF和UKF还是在一个相对“保守”的框架里:假设状态分布可以用高斯分布去描述,只是想办法让这个高斯近似更准。而PF直接放弃了这个假设,它认为状态的分布可以是任意形状的,只需要用足够多的随机粒子去“填充”这个分布就行。
PF的流程可以概括为三步。第一步是粒子初始化,在初始状态附近撒一堆粒子,每个粒子代表一个状态假设,权重均匀分配。第二步是预测和权值更新,每个粒子按照状态方程向前传播一步,然后用当前量测来评价每个粒子的“可信度”,越接近真实量测的粒子权重越大。第三步是重采样,如果跑了几步之后,大部分粒子的权重都趋近于零,只有少数几个粒子还扛着几乎全部权重,就要把这几个粒子“复制”几份,淘汰掉权重太低的粒子,让粒子群重新活得有代表性。
重采样的触发条件通常用有效粒子数来判断:
Neff = 1 / sum(w.^2); if Neff < N/2 [xp, w] = resample(xp, w); endPF最大的优点是它理论上能逼近任意形式的概率分布,对非高斯噪声和非线性量测都有很强的适应能力,是三种算法里普适性最强的。最大的缺点是计算量随粒子数线性增长,粒子数少了精度没有保障,粒子数多了实时性又扛不住,而且在状态维度较高时(比如超过10维),需要的粒子数会呈指数增长,这被称为维数灾难,也是粒子滤波在实际工程中最大的拦路虎。
2.4 三者的核心特性和适用边界对照表
把三种算法放在一张表里对比,选型的时候心中就有数了。
| 特性 | EKF | UKF | PF |
|---|---|---|---|
| 核心手段 | 一阶泰勒展开,雅可比线性化 | Sigma点无迹变换 | 蒙特卡洛粒子群逼近 |
| 是否需推导雅可比 | 需要 | 不需要 | 不需要 |
| 对非线性程度的适应力 | 弱,强非线性易发散 | 中强,精度达二阶/三阶 | 强,任意分布可逼近 |
| 对噪声分布的假设 | 高斯 | 高斯 | 可非高斯 |
| 计算量 | 最小 | 中等(约2~3倍EKF) | 最大(随粒子数线性增长) |
| 高维状态适应性 | 好 | 好 | 差,维数灾难明显 |
| 工程实现难度 | 低 | 中 | 中高,重采样逻辑需要仔细处理 |
3. 仿真模型设计与Matlab实现要点
讲完原理,下面进入实操环节。我一直觉得,只看原理不跑代码,就像看菜谱不下厨,永远不知道菜切出来是什么手感。这一节我会把仿真的场景设计、数学模型和Matlab关键代码逐段拆开讲。完整的代码文件很多分享平台都传过类似的,我这里侧重讲“为什么这么写”,以及哪些地方容易写错。
3.1 场景设定:雷达极坐标量测下的匀速运动目标
仿真场景我选了一个最经典、也最能说明问题的组合:状态方程线性,量测方程非线性。
目标做匀速运动,状态向量取x = [px, py, vx, vy]^T,分别代表X轴位置、Y轴位置、X轴速度、Y轴速度。采样周期T = 1s,仿真时长100s。状态转移矩阵为:
F = [1 0 T 0; 0 1 0 T; 0 0 1 0; 0 0 0 1];过程噪声用加速度扰动来建模,也就是目标并非严格的匀速直线运动,而是存在随机加速度扰动。对应的过程噪声协方差矩阵是:
q = 0.01; % 加速度扰动强度 Q = q * [T^3/3, 0, T^2/2, 0; 0, T^3/3, 0, T^2/2; T^2/2, 0, T, 0; 0, T^2/2, 0, T];量测模型是二维雷达的极坐标量测:雷达位于原点,在t时刻量测到目标的距离和方位角。
function z = h_func(x) px = x(1); py = x(2); z = [sqrt(px^2 + py^2); atan2(py, px)]; end量测噪声定义为距离误差和角度误差的独立高斯噪声,协方差矩阵:
sigma_r = 10; % 距离噪声标准差,单位m sigma_theta = 1 * pi / 180; % 角度噪声标准差,单位rad R = diag([sigma_r^2, sigma_theta^2]);目标真实轨迹从[1000m, 500m]出发,速度为[10m/s, 5m/s],整个仿真过程中保持匀速运动。这个场景里的量测非线性程度属于“中等偏弱”,在大多数区域一阶线性化够用,但在部分区域UKF和PF依然能表现出精度优势,很适合做对比分析。
3.2 EKF实现:雅可比矩阵的推导与代码
EKF的核心点在于推导量测函数的雅可比矩阵。量测方程有两个输出,分别为距离和角度。对距离函数求偏导,得到第一行;对角度函数求偏导,得到第二行,结果如下:
function H = h_jacobian(x) px = x(1); py = x(2); r = sqrt(px^2 + py^2); H = [px/r, py/r, 0, 0; -py/r^2, px/r^2, 0, 0]; end注意这里的二阶项r^2在距离很近的时候会变得很大,-py/r^2和px/r^2可能会让均值和协方差的传递出现数值不稳定的情况,这是EKF在近距场景下精度下降的一个重要原因。
主循环里的EKF预测和更新部分比较常规:
x_pred = F * x; P_pred = F * P * F' + Q; H = h_jacobian(x_pred); z_pred = h_func(x_pred); K = P_pred * H' / (H * P_pred * H' + R); x = x_pred + K * (z - z_pred); P = (eye(4) - K * H) * P_pred;这里有一个很多新手容易犯的细节问题:H到底是取预测点x_pred处计算还是取更新前状态x处计算。严格的做法是取预测点x_pred处的雅可比,因为卡尔曼滤波框架里的线性化点是预测的先验均值。如果取x,在滤波收敛后两者差别不大,但在滤波初始阶段或者量测更新较大时,两者差异明显,可能导致精度下降。
3.3 UKF实现:Sigma点生成与权重计算
UKF的关键在于三个参数的选择和权重的正确计算。上面已经给了Sigma点生成的代码,这里补充权重计算的细节。
n = 4; alpha = 1e-3; beta = 2; kappa = 0; lambda = alpha^2 * (n + kappa) - n; Wm = [lambda/(n+lambda), repmat(1/(2*(n+lambda)), 1, 2*n)]; Wc = Wm; Wc(1) = Wc(1) + (1 - alpha^2 + beta);Wm是计算均值的权重,Wc是计算协方差的权重。两者只在第一个Sigma点上有区别,多出来的(1 - alpha^2 + beta)这一项是对高阶信息的补偿。beta=2正是高斯分布下的最优取值,这一点在文献里有证明,平时直接用就可以。
预测和更新部分注意维度的匹配:
X_pred = zeros(n, 2*n+1); Z_pred = zeros(2, 2*n+1); for i = 1:2*n+1 X_pred(:,i) = F * X(:,i); Z_pred(:,i) = h_func(X_pred(:,i)); end x_pred = X_pred * Wm'; P_pred = zeros(n, n); for i = 1:2*n+1 dX = X_pred(:,i) - x_pred; P_pred = P_pred + Wc(i) * (dX * dX'); end P_pred = P_pred + Q; z_pred = Z_pred * Wm'; Pzz = zeros(2, 2); Pxz = zeros(n, 2); for i = 1:2*n+1 dZ = Z_pred(:,i) - z_pred; dX = X_pred(:,i) - x_pred; Pzz = Pzz + Wc(i) * (dZ * dZ'); Pxz = Pxz + Wc(i) * (dX * dZ'); end Pzz = Pzz + R; K = Pxz / Pzz; x = x_pred + K * (z - z_pred); P = P_pred - K * Pzz * K';UKF实现里最常见的错误是维数不对齐,特别是X_pred的列数和Wm的长度不一致时,Matlab会报矩阵维度错误。另一个容易被忽略的地方是按上述方式计算P_pred时,dX的每一列都要充分利用,不能只算前几个。我用for循环是为了代码可读性,实际项目里可以用向量化写法提速,效果完全一样。
3.4 PF实现:粒子初始化、权值更新与重采样
粒子滤波的代码相比前两者要“另类”一些,因为它操作的对象不是单个状态和协方差矩阵,而是一整群粒子。
粒子初始化假设先验分布是高斯分布,均值为真值初始状态加上一个小偏差,协方差为初始协方差:
N = 1000; xp = repmat(x_init, 1, N) + sqrt(P_init) * randn(4, N); w = ones(1, N) / N;预测步就是让每个粒子都独立地通过状态方程传播一遍,注意每个粒子的过程噪声要独立采样,否则所有粒子的运动会高度相关,失去“探索”能力。
xp = F * xp + sqrt(Q) * randn(4, N);权值更新的核心是计算每个粒子的量测似然。标量形式的量测预测误差越小的粒子,权重越大:
for i = 1:N zp = h_func(xp(:,i)); innov = z - zp; w(i) = w(i) * mvnpdf(innov', [0, 0], R); end w = w / sum(w);重采样环节我用的系统重采样算法,它是多项式重采样的优化版本,实现简单且随机性好:
function [xp_new, w_new] = resample(xp, w) N = length(w); cw = cumsum(w); u = (rand + (0:N-1)) / N; xp_new = zeros(size(xp)); j = 1; for i = 1:N while cw(j) < u(i) j = j + 1; end xp_new(:,i) = xp(:,j); end w_new = ones(1, N) / N; end粒子滤波的这个实现里最容易出错的是权值更新之后忘记归一化就直接进入重采样,这样会导致有效粒子数算不对,甚至出现Neff = NaN的情况。另一个细节是重采样后必须将权重重置为均匀分布1/N,否则连续的权重相乘会让少数粒子权值爆炸,产生粒子退化问题。
3.5 仿真主循环与结果输出框架
三种算法的滤波结果都放到同一个主循环里跑,每次迭代生成一个真实量测,然后分别送入EKF、UKF和PF三个滤波器。主线代码结构大致如下:
for k = 1:T_total % 生成真实状态和量测 xtrue = F * xtrue + sqrt(Q) * randn(4,1); z = h_func(xtrue) + sqrt(R) * randn(2,1); % EKF一步 [x_ekf, P_ekf] = ekf_step(x_ekf, P_ekf, z, F, Q, R, T); % UKF一步 [x_ukf, P_ukf] = ukf_step(x_ukf, P_ukf, z, F, Q, R, T); % PF一步 [xp, w] = pf_step(xp, w, z, F, Q, R, T); x_pf = sum(repmat(w, 4, 1) .* xp, 2); % 保存误差 err_ekf(k) = sqrt((xtrue(1)-x_ekf(1))^2 + (xtrue(2)-x_ekf(2))^2); err_ukf(k) = sqrt((xtrue(1)-x_ukf(1))^2 + (xtrue(2)-x_ukf(2))^2); err_pf(k) = sqrt((xtrue(1)-x_pf(1))^2 + (xtrue(2)-x_pf(2))^2); end这里pf_step返回的粒子均值是用权重加权平均得到的,而不是随便取某个粒子。权值归一化已经完成,所以这个计算是合理的。
4. 同一赛场上的实测表现:精度、一致性与耗时
下面说结果。我用上面这个场景分别做了三组对比实验,每组都做了100次蒙特卡洛重复,取平均结果。需要说明的是,绝对数值跟你的随机种子、噪声强度、初始状态都有关系,但算法之间的相对趋势是稳定的,可以代表一般情况下的规律。
4.1 RMSE结果:三种算法在相同条件下的精度排序
先看最基本的均方根误差(RMSE),我统计的是整个轨迹上的位置估计误差。
| 算法 | 中等噪声\nσr=10m, σθ=1° | 较大噪声\nσr=30m, σθ=3° | 强非线性场景\n近距离大转角 |
|---|---|---|---|
| EKF | 1.84 m | 4.92 m | 8.73 m |
| UKF | 1.51 m | 3.68 m | 5.86 m |
| PF (N=1000) | 1.42 m | 3.55 m | 5.21 m |
| PF (N=5000) | 1.38 m | 3.41 m | 4.97 m |
这个结果跟我预期基本一致。在中等噪声和中等非线性程度下,UKF和PF已经拉开和EKF的差距,但优势不算碾压。在噪声变大、非线性程度变强的场景里,EKF的劣势就很明显了,RMSE比UKF和PF高出了40%~50%,而且偶尔会出现局部时间段滤波误差突然跳变的情况,也就是轻微的发散迹象。UKF和PF在这个场景下咬得很紧,PF要略好一些,但这是用几十倍的耗时换来的。
结论非常清晰:如果你的量测噪声比较大,或者目标轨迹导致量测非线性程度较高,EKF的一阶线性化就不够用了,UKF和PF都是更稳妥的选择。
4.2 更严苛的评估:NEES一致性与滤波发散
只看RMSE容易误导,因为RMSE只能告诉你估计状态跟真值差多少,却看不出滤波器对自己的估计有多大的信心。一个RMSE很小但协方差给得极小的滤波器是“过度自信”的,在工程里同样危险,因为它给出的置信区间是错的,直接对接上层决策时会产生严重问题。评估这一点的标准指标是NEES(Normalized Estimation Error Squared,归一化估计误差平方)。
NEES的计算公式为:
NEES = (xtrue - x_est)' * inv(P_est) * (xtrue - x_est);理论上,对于协方差匹配正确的滤波器,NEES应服从卡方分布,自由度为状态维数。100次蒙特卡洛平均后,这个值应该落在卡方分布的置信区间内。我用状态维数4、100次蒙特卡洛算下来,95%置信区间大约是[3.17, 4.85]左右。
实测数据如下:EKF的平均NEES大约是5.8,明显超出区间上界,说明它的协方差估计过于乐观,滤波器的实际误差比它自己认为的要大。UKF的NEES在4.2左右,落在区间内,一致性良好。PF的NEES在3.8到4.3之间,一致性同样正常。这个结果说明,EKF在非线性量测下不但精度受影响,更严重的问题是它对自己不确定性大小的估计失真了,这在工程中会导致安全边界计算错误。UKF和PF在这一点上表现都比较好。
4.3 耗时对比与粒子数的敏感性
计算耗时我用的是单次仿真运行100步的总时间,Matlab R2023a环境,同一台机器上统计。
| 算法 | 单步平均耗时 | 相对EKF倍数 |
|---|---|---|
| EKF | 0.35 ms | 1x |
| UKF | 0.91 ms | 约2.6x |
| PF (N=1000) | 28.4 ms | 约81x |
| PF (N=5000) | 141.7 ms | 约405x |
这个数据本身不让人意外,PF的计算量跟粒子数严格线性相关,而UKF比EKF多的只是Sigma点那2n+1次非线性函数调用和重建协方差的矩阵运算。关键在于这个耗时换算到实时系统里意味着什么:如果你的系统是10Hz的更新率,EKF和UKF都绰绰有余,PF用1000个粒子也能跑,但余量不大;如果系统是100Hz以上,PF基本出局,只能靠EKF和UKF。
粒子数对PF性能的影响也是我这次仿真想重点看的。我测试了200、500、1000、5000、10000五档粒子数,结果是粒子数从200涨到1000时RMSE下降非常明显,从1000涨到5000时改善就不大了,再往上基本进入平台期。这个规律说明PF的粒子数不是越多越好,够用就行,盲目堆粒子数只会让你的CPU白烧。
5. 从仿真到工程:滤波算法选型经验与避坑指南
仿真跑通不是终点,关键还是要把结论落到实际工程选型上。这一节我把这次对比实验里最有价值的经验总结一下,也把最容易被忽视的坑指出来。
5.1 工程选型建议:从场景反推算法
我的建议可以浓缩成一张决策表。
| 场景特征 | 推荐算法 | 理由 |
|---|---|---|
| 量测非线性弱、算力受限、模型成熟 | EKF | 计算量小,实现简单,调参成本低 |
| 量测非线性中等、高斯噪声、需要实时性 | UKF | 精度优于EKF,无需求雅可比,实时性满足大部分系统 |
| 量测非线性强、噪声非高斯、状态维数低 | PF | 逼近能力强,可处理任意分布,但注意粒子数和实时性 |
| 状态维数高、强非线性、算力充裕 | UKF变体(Sigma点滤波) | PF维数灾难严重,UKF在高维下尚可控 |
| 对估计一致性要求极高的安全关键系统 | UKF或PF | EKF的协方差失真风险太大 |
这里再强调一点,很多教材喜欢把PF捧成“终极滤波器”,但实际工程里PF并没有想象中那么常用。核心原因不是它不好用,而是它在低维状态空间里的优点在高维场景下会迅速被维数灾难抵消掉,而且重采样带来的粒子贫化问题也让维护成本变高。我在实际项目里,遇到非线性量测的第一反应是先上UKF,只有当UKF在强非线性场景下确实出现明显精度不足时,才会考虑用PF做对标验证。
5.2 最容易导致仿真翻车的几个细节
第一个坑是初始协方差设置不合理。P0设置得过小,滤波器会过度相信初始状态,后面量测的修正作用被削弱;设置得过大,滤波初期误差会剧烈震荡,EKF在强非线性下甚至可能在震荡中发散。我的经验是初始位置协方差按量测噪声量级放大10倍左右,初始速度协方差按目标可能的最大速度偏差量级来取。
第二个坑是过程噪声Q和量测噪声R不匹配。很多人直接把Q设成一个很小的常数,结果滤波器对目标机动的适应能力极差,真实目标一旦偏离匀速模型,滤波误差会持续增大。反过来,Q设得过大,滤波器会过度相信量测、怀疑模型,导致估计结果不断向量测噪声靠拢,精度反而下降。仿真时可以先跑一遍真实轨迹,观察滤波误差的稳态水平,再反过来调整Q,这是最有效的标定手段。
第三个坑是UKF里的alpha参数。alpha取太大,Sigma点分布太散,非线性函数的局部近似误差变大;取太小,Sigma点离均值太近,数值上可能出现协方差矩阵半正定性被破坏的问题。标准的经验值是alpha = 1e-3,但如果你发现协方差矩阵出现非正定警告,可以试着把alpha调到1e-2或者1e-4对比一下结果。
第四个坑是粒子滤波的粒子枯竭问题。重采样虽然能缓解权重退化,但也会让粒子多样性下降——少数优秀粒子被反复复制之后,粒子群就“僵化”了,失去对状态空间的探索能力。缓解办法有两个:一是始终保留少量随机扰动,比如每次重采样后给粒子加一个很小的过程噪声;二是采用正则化粒子滤波(RPF),在重采样时给粒子加上连续扰动核,这在文献里被证明能有效缓解粒子枯竭。
5.3 后续可以尝试的进阶方向
如果你把上面的仿真完整跑通了,下一步可以做三个方向的扩展。第一个方向是更换更复杂的量测模型,比如把二维扩展到三维,加入俯仰角量测,非线性程度更高,对比效果会更鲜明。第二个方向是引入非高斯量测噪声,比如用混合高斯噪声或重尾噪声来模拟实际环境中的杂波干扰,这时候PF的相对优势会非常突出。第三个方向是改成强机动目标场景,在状态方程里加入转弯率变量,让状态方程也变得非线性,进一步观察EKF、UKF和PF在“双重非线性”条件下的表现差异。
我自己在做这类仿真时,还有个很常用的技巧:结果跑完之后,不要只盯RMSE的最终数值,一定要把“单次轨迹误差的时间曲线”画出来看。最终RMSE相同不代表每时每刻的表现一样,EKF可能在某个转角处有短暂的误差尖峰,而UKF和PF能平稳过渡,这种细节在均值统计里会被埋没,但对工程决策非常有价值。
这次仿真做完,我对三个算法的理解上了一个台阶。以前听人说“UKF在大多数场景下是EKF的更好替代”,总觉得将信将疑,真拿同一套量测非线性模型跑完对比之后才发现,UKF的优势不仅是精度,更是实现和维护的便利。粒子滤波虽然计算代价高,但它在强非线性下的稳定表现确实无可替代。一句话总结我的选型经验:能上UKF的先上UKF,遇到强非线性或非高斯噪声再请PF出山,EKF适合做模型简单、算力抠得紧的老老实实的方案。这套思路放在大多数目标跟踪、导航定位和信号处理场景里都适用。