最近把手头这个弹道目标跟踪仿真项目重新整理了一遍，感觉有很多值得记录的东西。这个项目的任务很明确：对一个受到空气阻力影响的弹道目标进行状态估计，状态量包含高度、速度、弹道系数，分别用扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）两条技术路线去实现，全程在Matlab环境下仿真验证，并附完整代码。弹道目标的状态估计在雷达数据处理、飞行器轨迹预测、再入段跟踪这些方向属于非常经典的课题，很多实验室和研究团队都会拿它当卡尔曼滤波系列的练手项目。但"入门级"三个字并不等于没坑可踩，真正从公式走到能跑的代码，中间涉及建模假设、坐标约定、雅可比推导、数值稳定处理等一系列细节。这篇博文就按我这个项目的实际推进顺序来写：先说清楚弹道目标是怎么建模的，再分别拆EKF和UKF的实现要点，然后给出两者在同一仿真场景下的对比结果，最后把代码结构和调试经验一并交代，希望能帮你少走几个弯路。

1. 弹道目标的运动建模：空气阻力不是装饰品

1.1 为什么状态向量必须包含弹道系数

很多刚接触目标跟踪的读者会有一个习惯性思维：目标运动无非就是匀速直线运动模型或者匀加速运动模型，状态量取位置和速度就够了。这个思路在近距匀速运动目标上没问题，但遇到弹道目标就会露馅。

弹道目标飞行过程中受重力和空气阻力两个主要力的作用。重力加速度基本恒定（在几十公里高度范围内变化不大），而空气阻力的大小和空气密度、速度平方、物体本身的气动特性强相关。同样的速度下，一个质量大、阻力小的弹丸和一个质量小、迎风面积大的物体，轨迹衰减速度完全不同。所以只估计高度和速度，模型里就没有足够的自由度去描述"这个目标到底有多容易被空气减速"这一物理特性，滤波结果很快就会偏离真值。

工程上习惯用弹道系数β来综合描述目标的空气动力学属性：

β = m / (Cd * A)

其中m是质量，Cd是阻力系数，A是参考迎风面积。β的量纲是kg/m²，它直观的含义是"单位迎风面积上分配了多少质量"。β越大，说明目标越"重"、越不容易被空气阻力减速；β越小，说明目标"轻飘飘"、阻力影响显著。

所以状态向量取x = [h, v, β]^T不是随便拍脑袋定的，它对应的是"目标当前在哪、动得多快、气动特性如何"三组信息，三者之间有内在的动力学耦合关系。只测高度和速度而不估计β，模型就无法区分"速度衰减是因为高度越高空气越稀薄"还是"这个目标本身阻力就大"。

1.2 从连续运动方程到离散状态方程

垂直方向弹道运动的连续时间模型可以写成：

dh/dt = v

dv/dt = -g - ρ(h) * v * |v| / (2β)

dβ/dt = 0

速度方程里为什么要用v * |v|而不是v²？因为v的符号是有物理意义的。我定义的v以向上为正，上升段速度向下取负值，下降段速度向上为正（如果目标是下落体）时，阻力方向要始终与速度方向相反。用v*|v|这个形式，阻力项自动保持"始终阻挡运动"的特性，比直接写v²更严谨。若目标全程是上升段，v>0恒成立，两者等价；但要做完整的再入段或抛射体轨迹仿真时，sign问题就无法回避了。

空气密度ρ(h)我采用指数大气模型：

ρ(h) = ρ₀ * exp(-h/H)

ρ₀取1.225 kg/m³，H取8500m。这个模型的含义是：高度每上升约8.5公里，空气密度下降一个自然指数倍。它比常数密度模型复杂一个量级，但完全是值得的——高空稀薄空气对弹道目标阻力的影响极其显著，如果把ρ当常数，滤波结果在高度跨度大的场景下基本不可用。

连续模型不能直接用于计算机离散仿真和卡尔曼递推，需要离散化。最朴素的做法是欧拉法，一步到位：

h_{k+1} = h_k + v_k * dt

v_{k+1} = v_k + [-g - ρ(h_k) * v_k * |v_k| / (2β_k)] * dt

β_{k+1} = β_k

这里有个重要问题：生成真值轨迹时如果也用欧拉法，等于"用模型生成数据再用模型滤波"，属于循环论证，无法暴露模型离散误差。我在项目中生成真值用的是Matlab的ode45求解器，设置相对误差1e-8，把连续方程积分到各采样时刻，得到高精度参考轨迹。这样滤波过程中使用的离散状态方程只是"近似描述"，真值和模型之间的失配会更接近实际工程情况。

状态方程的非线性源有两个：一个是v*|v|项带来的二次非线性，另一个是ρ(h)对高度的指数依赖。这两个非线性源恰恰是EKF和UKF要处理的难点，也是后面分析两者差异的主战场。

观测模型方面，我用了雷达对高度和速度的直接观测，观测方程是线性的：

z_k = H * x_k + v_k

H = [1 0 0; 0 1 0]

观测噪声v_k服从零均值高斯分布，标准差分别设为σ_h和σ_v。这里观测是线性的，过程是非线性的，属于状态估计里很典型的一种混合结构。实际雷达一般输出距离和角度，需要经过坐标转换才能得到高度和速度的观测量，但坐标转换不属于本次滤波算法的研究重点，所以做了简化处理。如果你要直接套用真实雷达数据，需要在观测方程里额外考虑坐标转换的雅可比矩阵。

2. EKF在弹道跟踪中的实现细节

2.1 雅可比矩阵推导的完整过程

EKF的核心思想是"线性化"：在每次预测时，把非线性状态方程在当前估计值附近做一阶泰勒展开，用雅可比矩阵代替线性系统中的状态转移矩阵。思路本身不复杂，但具体到弹道目标这个模型，雅可比矩阵每一个元素的推导都需要仔细校验，算错一个符号滤波就会发散。

设状态向量的三个分量f₁、f₂、f₃分别是离散状态方程的右端：

f₁ = h + v*dt

f₂ = v + [-g - ρ(h) * v * |v| / (2β)] * dt

f₃ = β

状态转移雅可比F的3×3个元素按以下方式计算。

第一行：∂f₁/∂h = 1，∂f₁/∂v = dt，∂f₁/∂β = 0。

第二行是重头戏。∂f₂/∂h涉及密度对高度的导数：

∂f₂/∂h = -[∂ρ/∂h * v * |v| / (2β)] * dt

由于ρ(h) = ρ₀ * exp(-h/H)，∂ρ/∂h = -ρ(h)/H，代入后：

∂f₂/∂h = dt * ρ(h) * v * |v| / (2β * H)

这里密度项对高度的负导数被消掉了，最终得到一个正值——高度越高，空气密度越小，阻力对速度的影响越小，所以高度对速度增量的偏导是正的。初学者很容易在这里落下负号，导致雅可比矩阵物理意义错误。

∂f₂/∂v = 1 - dt * ρ(h) * |v| / β

这里用到的关键公式是d(v*|v|)/dv = 2|v|。阻力项对v的导数是负的，代表速度越大阻力越大、速度增量越小的负反馈机制。

∂f₂/∂β = dt * ρ(h) * v * |v| / (2β²)

这里分母出现β²，来自对1/β求导。这个元素的意义是：弹道系数越大，阻力越小，速度衰减越弱，所以β对速度增量的偏导为正。

第三行全零，除了∂f₃/∂β = 1。

得到雅可比矩阵后，EKF的预测和更新就按标准流程走：

% EKF预测 x_pred = zeros(3,1); x_pred(1) = x(1) + x(2)*dt; x_pred(2) = x(2) + (-g - rho(x(1))*x(2)*abs(x(2))/(2*x(3)))*dt; x_pred(3) = x(3); F = ekf_jacobian(x, dt, g, rho0, H_scale); P_pred = F * P * F' + Q; % EKF更新 z_pred = H * x_pred; S = H * P_pred * H' + R; K = P_pred * H' / S; x = x_pred + K * (z - z_pred); P = (eye(3) - K * H) * P_pred; P = (P + P') / 2;

2.2 滤波初始化和协方差矩阵的工程处理

EKF的初始状态x₀和初始协方差P₀怎么取，直接决定滤波能否收敛。我在仿真里是这样设置初始值的：真实状态是h₀、v₀、β₀，滤波初始值取真实值加一个随机偏差，初始协方差P₀根据这个偏差的大小来定。

这里要解释一个初学者容易踩的坑：P₀不是随便给个对角阵就完事的。P₀的元素反映的是你对初始状态的不确定程度，如果初始偏差为Δh、Δv、Δβ，那么P₀对角元大致取(3Δ)^T级别的量级。取太小会让滤波器"过度自信"，后续观测来不及修正，滤波结果被冻结在初始偏差附近；取太大则前期增益过大，滤波初期容易震荡甚至发散。

还有一个细节：速度初值v₀与弹道系数初值β₀之间存在很强的耦合。初始速度差一点会导致阻力变化，阻力变化又会影响速度衰减估计，因此设计P₀时不要把所有元素当成相互独立的。更稳健的做法是给P₀的非对角元赋一个小的交叉项，或者在过程噪声Q里给速度分量留出足够的余地，让滤波器自己调整。

另一个工程问题是P矩阵失去对称正定性。EKF递推过程中因为数值舍入误差，P可能逐渐变得不对称，甚至出现负特征值，导致后续Cholesky分解直接报错。我在代码里做了一步保护：每次更新完P后强制对称化，并对角线加一个极小量：

P = (P + P') / 2; P = P + 1e-12 * eye(3);

这个处理看起来不起眼，却能挡住90%以上"滤波器跑着跑着就崩了"的问题。UKF里这个保护更加重要，因为sigma点生成依赖协方差矩阵的平方根分解。

3. UKF的无迹变换策略与代码落地

3.1 sigma点参数选择的门道

UKF避开了解析求导，改用无迹变换的思想：在状态分布中挑选一组确定性采样点（sigma点），让这些点经过非线性函数传播后，用加权统计估计均值和协方差。这个思路对非线性强度的容忍度比一阶线性化高得多，尤其适合弹道目标这种包含速度平方和指数密度项的问题。

对三维状态，取2n+1 = 7个sigma点。核心参数有三个：

• α：决定sigma点相对均值的散布程度，典型范围1e-3到1
• β：用于融入状态分布的先验信息，高斯分布时最优取2
• κ：辅助缩放因子，通常取0，或满足n+κ=3

缩放参数λ = α²(n+κ) - n，gamma = sqrt(n+λ)。sigma点列生成如下：

lambda = alpha^2 * (n + kappa) - n; gamma = sqrt(n + lambda); X = zeros(n, 2*n+1); X(:,1) = x; xP = gamma * chol(P)'; % 平方根分解 for i = 1:n X(:, i+1) = x + xP(:, i); X(:, i+n+1) = x - xP(:, i); end

权重分配：

Wm(1) = lambda / (n + lambda); Wc(1) = lambda / (n + lambda) + (1 - alpha^2 + beta); for i = 2:2*n+1 Wm(i) = 1 / (2 * (n + lambda)); Wc(i) = 1 / (2 * (n + lambda)); end

实际调参时，α的选择最影响效果。把α取到接近1，sigma点散布范围大，对强非线性系统的鲁棒性更好；但代价是采样点离均值太远，当状态分布比较窄时可能采样到概率极低的区域。α取1e-3时sigma点集中在均值附近，对弱非线性更精准，但协方差平方根计算时数值精度要求更高。我最终在弹道模型上用的α=0.6，兼顾两者。当然这只是本项目的经验值，换场景得重新调。

3.2 UKF实现中容易写错的几个地方

UKF相比EKF公式不算难，但有几个细节问题我折腾了不少时间。

第一，所有sigma点必须通过同一个非线性状态方程传播。这个说法听起来像废话，但实际代码里很容易因为向量化处理失误，导致不同sigma点用了不同的参数。比如ρ(h)的计算，每个sigma点的高度分量不同，密度自然不同，必须逐点计算。初学者可能会图省事用均值处的密度代替所有sigma点的密度，这会破坏sigma点的一致性，滤波精度大打折扣。

第二，观测更新的协方差计算中，观测预测值z_pred的均值要由所有sigma点的观测预测加权得到，而不是用x_pred对应的观测。这两者的差异虽然微小，但在状态和观测强相关时会明显影响滤波增益的优化。

第三，处理观测更新时同样存在正定性问题。过程更新得到的P_pred经过观测更新后可能出现非正定，我的做法是每次状态更新和协方差更新完成后都做对称化处理。

UKF完整更新循环的核心片段：

% 过程传播 Y = zeros(n, 2*n+1); for i = 1:2*n+1 Y(1,i) = X(1,i) + X(2,i)*dt; Y(2,i) = X(2,i) + (-g - rho(X(1,i))*X(2,i)*abs(X(2,i))/(2*X(3,i)))*dt; Y(3,i) = X(3,i); end x_pred = sum(repmat(Wm, n, 1) .* Y, 2); P_pred = zeros(n, n); for i = 1:2*n+1 dx = Y(:,i) - x_pred; P_pred = P_pred + Wc(i) * (dx * dx'); end P_pred = P_pred + Q; % 观测更新 Z_pred = H * Y; z_pred = sum(repmat(Wm, m, 1) .* Z_pred, 2); Pzz = zeros(m, m); Pxz = zeros(n, m); for i = 1:2*n+1 dz = Z_pred(:,i) - z_pred; dx = Y(:,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';

注意到我用H * Y，其中Y是过程传播后的sigma点矩阵。这一步的整体观测量是由每个sigma点的观测量加权综合得到的，和线性观测模型下先算均值再观测在数学上等价，但代码实现时用前者更通用——如果未来把观测模型换成非线性的测距/测角方程，这套代码可以直接复用，不需要改结构。

4. 同一弹道场景下EKF vs UKF：仿真对比怎么看

4.1 实验设置：蒙特卡洛与评价指标

算法做好之后，单次仿真结果随机性太强，不能作为评判依据。我在项目里设置了这样一个基准场景：

• 目标初始高度h₀ = 30km，初始速度v₀ = 1800m/s（上升段）
• 弹道系数β = 6000 kg/m²
• 采样周期dt = 0.5s，仿真时长60s
• 观测噪声：σ_h = 100m，σ_v = 20m/s
• 初始估计偏差：高度偏差200m，速度偏差60m/s，弹道系数偏差20%（初始β估计为4800）

这套参数模拟的是一个高空高速目标，空气密度随高度迅速变化，过程模型非线性较强。然后跑了200次蒙特卡洛仿真，每次观测噪声和初始估计偏差随机生成，统计两种算法的均方根误差（RMSE）和卡尔曼滤波收敛情况。

RMSE的计算方式：

rmse = sqrt(mean((x_est - x_true).^2, 2));

分别对高度、速度、弹道系数三个状态量计算。同时记录单步平均耗时，用Matlab的tic/toc统计。

4.2 对比结果的可视化解读

先说结论：UKF在精度上全面占优，尤其是在弹道系数β的估计上优势明显。EKF虽然也能收敛，但在非线性较强的前期阶段出现了明显的偏差振荡。

高度估计的RMSE两者差别不算大，EKF大约比UKF高出15%到20%，都在可接受范围内。这是因为高度观测存在且观测噪声较小，状态量的可观测性较强，即使过程模型线性化误差存在，观测也能把估计拉回真值附近。

速度估计开始拉开差距。这是因为速度是通过运动方程间接估计的，而速度方程里的非线性项v*|v|对雅可比线性化非常敏感。EKF的一阶近似在速度方向上不够精确，表现为滤波初期速度估计出现系统偏差，直到观测持续修正后才慢慢消除。UKF的sigma点传播保留了高阶统计信息，速度方向上的偏差小得多。

差距最明显的是弹道系数β。EKF的β估计曲线在最初10秒内出现了一个明显波动，有时候甚至偏离20%左右才被拉回来；UKF的β估计则相对平滑，约5秒后就稳定在真值附近。这个结果在物理上可以解释：β是通过速度衰减率间接观测的，速度衰减率本身是阻力的函数，阻力又同时依赖密度、速度和β三个量。要从中解耦出β，对过程模型的精度要求极高。EKF的一阶截断误差在这个间接观测环节被放大，导致β的"反应速度"变慢。

计算耗时方面，UKF单步要传播7个sigma点，每一轮时间大概是EKF的3到4倍。在我的测试机器上EKF单步约0.3ms，UKF约1.1ms。不过两者都远小于0.5s的采样周期，所以实时性都不构成瓶颈。

如果只做一次仿真，随机性可能掩盖真实差距；但统计学上200次蒙特卡洛的平均结果趋势是很明确的。这个结论也提醒你：EKF并没有到"过时"的程度，如果系统非线性较弱、模型匹配较好，EKF的性价比很高；但像弹道目标这种强非线性场景，UKF的额外计算成本是值得花的。

5. Matlab仿真代码的结构设计与调试经验

5.1 代码模块划分与关键函数

完整代码包我按功能拆成了五个文件，这样的好处是算法、场景、后处理完全解耦，想换一个观测模型或者改仿真场景不用动核心滤波函数。

主脚本的核心流程就三段：先生成轨迹和观测，然后分别调用EKF和UKF，最后统一绘图对比。我特意把两个滤波函数设计成完全相同的输入输出格式，这样就可以在同一套后处理代码里比较它们，也方便以后替换成其他滤波器比如粒子滤波。

弹道系数初始值对滤波收敛速度影响很大，主脚本里我把它设计成可以外部传入的参数，方便做敏感性分析。在实际运行中，你可能想测试"初始弹道系数偏差50%会怎样"，这时只需要改一个参数重跑主脚本，不用动滤波核心代码。

5.2 参数调节和常见坑

这个项目的调试过程其实花了大半的时间在参数和细节上。如果你准备跑通并改写这套代码，以下几条经验应该能帮你省点时间。

过程噪声Q的选取是整个滤波器能否长期稳定运行的关键。我的做法是让Q的对角元与状态量范围成比例，同时给β通道留一个非常小的过程噪声，代表"弹道系数在实际飞行中不完全恒定"的扰动。β的过程噪声取太大会让协方差无限膨胀，导致滤波器"忘掉"之前累积的信息，β估计出现随机游走；取太小则滤波器对β的变化不敏感，飞行器机动或质量变化时无法追踪。我的折中方案是让β的Q值约为初始β平方的1e-6量级。

观测噪声协方差R相对好取，直接按传感器标称精度填写即可。但注意如果你把观测从"高度+速度"改成"只测高度"来验证算法极限，R的取值要重新审视。只测高度时，速度v和β成为间接估计量，可观测性下降，滤波容易发散。遇到这种情况不要急着怀疑算法，先检查R是否合理，以及P₀中速度和β的不确定性是否过小。

UKF的Cholesky分解报错是最高频的故障。几乎所有"posdef error in chol"都源于P矩阵不正定。我的排查顺序是：先看P是否对称，再看P的特征值是否有负值，最后检查Q是否太小导致预测协方差奇异。对称化处理要在每一步做，不是只在初始化时做。

最后还有一个容易被忽略的细节：空气密度模型用的是指数大气，当高度降到接近海平面时ρ趋近ρ₀，但当高度为负（比如目标落地点低于参考面）时ρ会异常增大。我在密度函数里加了一个高度下限保护，h小于0时直接取ρ₀。这个保护在正常弹道场景下不会被触发，但在蒙特卡洛仿真中某些发散轨迹可能产生负高度，没有保护的话算法会直接崩溃。

整套代码的完整m文件和参数配置文件我打包在代码包中，里面有详尽的注释和运行说明，直接用Matlab打开运行main_ballistic_ekf_ukf.m就能复现出本文图中的对比曲线。

6. EKF和UKF代码对比后的总结与选型建议

写到这里，我想把自己做这个项目过程中体会最深的几点沉淀下来。回看整个建模和调试流程，弹道目标状态估计的核心难点其实不在滤波算法本身，而在"模型是否真实反映了物理过程"和"滤波器的假设是否与模型失配程度兼容"这两个层面。EKF和UKF在同一场景下的表现差异，根源也恰恰在这里。

从适应场景的角度说，如果你的系统模型非线性较弱、初始误差不大，EKF已经足够可靠，它的代码简单、计算量小，而且因为有明确的雅可比矩阵，调试时能更直观地定位问题。UKF虽然精度高，但sigma点采样的统计特性、参数调节的复杂度都会引入新的可变因素，在强非线性场景下收益才明显。像弹道目标这种带指数密度依赖、速度平方项、且需要从间接观测量中估计弹道系数的系统，UKF是更稳妥的起点。

如果你想把这套代码扩展到二维弹道跟踪，比如状态量增加水平距离x和水平速度vx，关键的两处是：状态方程里要增加水平方向的运动方程，观测模型要改成雷达距离-方位角格式，对应的雅可比矩阵和sigma点数都要跟着变。过程模型变成五维后sigma点从7个变成11个，计算量上升但UKF代码结构不需要改。这个扩展方向我建议你做完基础实验后去尝试一下，对理解状态估计的通用性会很有帮助。

如果准备在这个方向上继续深入，后面还可以考虑自适应噪声协方差、交互式多模型（IMM）对机动弹道目标的切换滤波，以及把仿真数据和真实雷达记录做对比验证。滤波算法的价值最终还是要落在真实系统的稳定运行上，仿真只是第一步。我个人的体会是，在仿真里把每一个参数的物理意义想清楚，远比盲目调参得到好看的曲线更有价值，因为前者给了你对付新问题的能力，后者只给了你一篇能用的报告。