基于EKF/UKF的路面附着系数估计：Matlab/Simulink仿真与调试

去年某底盘控制项目里，我们在一段压实雪路面上反复标定ABS逻辑，同一套滑移率门限在干沥青和冰雪上完全是两套结论。问题根子很明确：路面附着系数永远没法用传感器直接测。做路面附着系数估计，本质上就是用车轮、车身这些能测到的信号，反向推算轮胎和路面之间到底还剩多少抓地力。这篇我把自己在 Matlab/Simulink 里实现 UKF 和 EKF 做路面附着系数估计的全过程捋一遍，包括模型怎么建、滤波器怎么改、仿真怎么搭、参数怎么调，以及几个容易让人卡住半天的坑。如果你正准备做类似的车辆状态估计课题或底盘控制算法预研，这篇可以直接当成起步参考。

1. 为什么要在在线估计路面附着系数：ABS测试里的现实教训

1.1 路面摩擦的“瞬变”属性与固定参数的死穴

附着系数是轮胎-路面接触面能传递的最大切向力与法向载荷的比值，工程里常用 μ 表示。干燥沥青路面上 μ 峰值能到 0.8-0.9，压实雪路面立刻掉到 0.25 以下，冰面更是只有 0.1 上下。更麻烦的是，哪怕同一条路，刚下过雨、路面有积水、轮胎温度不同，μ 都会变。底盘控制里 ABS、TCS、ESP 全部依赖这个数来决策：到底允许车轮产生多大滑移率，或者应该给制动轮缸施加多大的压力。

早些年很多 ABS 方案采用固定门限加分段修正，标定时反复调整滑移率阈值，目的就是让控制逻辑在各路面条件下都能正常工作。问题在于固定参数是为某一类路面优的，换路面后要么过早介入、要么介入太晚。在我们那次冰雪测试里，干沥青标定出来的参数上雪地之后，ABS 介入明显滞后，某测试车差点冲出弯道。事后复盘大家都同意一件事：如果控制算法能实时知道当前路面附着系数，很多标定死角可以直接化解。

1.2 为什么不直接查表或用力相除

最直觉的做法是用纵向力除以法向载荷：μ = Fx / Fz。Fx 可以通过驱动扭矩或制动压力换算，Fz 大致等于车重分配。这个方法听上去简单，做起来全是噪声。Fx 的估计本身就有误差，制动压力到夹紧力的响应还有时间常数，纵向力没有充分激励的时候除出来一个毫无意义的数。还有更麻烦的：车轮在纯滚动状态时，轮胎纵向力很小，这时候把力除以载荷得到的根本不是附着上限，只能表示“当前用到的力”，离极限还差得远。

所以工程上更推荐用状态观测思路：把 μ 当成系统里的一个未知慢变量，和车速、轮速一起放进状态向量，用卡尔曼类滤波器同时估计。核心逻辑是，车轮动力学方程把纵向力 Fx 和 μ 联系起来，而 μ 的变化比车速慢得多，可以用一个随机游走模型近似。这就是 EKF/UKF 能干活的原因：它不要求你直接测 μ，只需要把可测信号在时间上的变化规律和车辆模型的物理关系揉进贝叶斯更新框架里。

2. 单轮车辆模型与 Burckhardt 轮胎模型：把物理过程写进状态方程

2.1 整车到单轮的降维：纵向控制问题没必要上十几自由度模型

路面附着系数估计问题里，如果只关心车辆纵向动力学（制动或驱动工况），用单轮模型就够了。十几自由度的动力学模型当然能描述整车姿态，但对估计器来说意义不大，还增加计算负担，不利于后续嵌入实时控制器。

单轮模型做以下假设：车辆直线行驶或侧偏角很小，忽略横向载荷转移；四个车轮受力近似对称，用一个等效车轮代表全车；轮胎垂直载荷近似不变；阻力简化为可以忽略或并入过程噪声。在这些假设下，两个方程构成系统骨架：

车轮旋转动力学方程：

I_w · dω/dt = T_d - T_b - F_x · R

整车纵向动力学方程：

m · dv/dt = F_x

其中，I_w 是车轮转动惯量，ω 是车轮角速度，R 是轮胎滚动半径，v 是车辆纵向速度，F_x 是轮胎纵向力，T_d 和 T_b 分别是驱动力矩和制动力矩。

纵向力由附着系数和垂直载荷构成：

F_x = μ(λ) · F_z

F_z 取单轮垂直载荷，μ 是滑移率 λ 的函数。制动工况下纵向滑移率定义为：

λ = (v - ω·R) / v

这里要特别提醒：滑移率分母用的是车辆速度 v，不是轮速。λ 的范围是从 0（纯滚动）到 1（车轮完全抱死，ω 为 0）。

2.2 Burckhardt 轮胎模型与“慢变参数”假设

纵向力模型我推荐 Burckhardt 模型，原因是它参数少、曲线形状可以直接反映不同路面的附着特性，而且后面做 EKF 求导时表达式也相对好处理。公式是：

μ(λ) = C1 · (1 - exp(-C2·λ)) - C3·λ

不同路面的 C1、C2、C3 系数差异明显：

• 干沥青：C1 ≈ 1.28，C2 ≈ 23.99，C3 ≈ 0.52
• 湿沥青：C1 ≈ 0.86，C2 ≈ 33.82，C3 ≈ 0.35
• 压实雪：C1 ≈ 0.19，C2 ≈ 94.13，C3 ≈ 0.06

这个模型的含义用大白话讲：随着滑移率增大，μ 先是快速上升，到达峰值后缓慢下降。峰值位置代表轮胎纵向力的极限，ABS 工作区间通常就在峰值附近。

现在问题来了：实际路面属于哪一类系数组合，估计器不知道。如果把 C1、C2、C3 全部作为待估计状态，系统维数上升且参数耦合严重，很难收敛。工程上的做法是把 μ 本身直接作为状态量，假设它在相邻采样时刻基本不变，用随机游走模型描述其动态：

dμ/dt = 0 + w_μ

w_μ 是过程噪声。这样虽然丢失了 μ 与滑移率之间的精确映射关系，但保留了估计器最需要的核心信息：当前路面的附着能力。车辆状态向量选为：

x = [v, ω, μ]ᵀ

整套系统的离散状态转移可以近似写为：

• v(k+1) ≈ v(k) + Δt · (F_x / m)
• ω(k+1) ≈ ω(k) + Δt · (T_d - T_b - F_x·R) / I_w
• μ(k+1) ≈ μ(k)

注意，我这里没有把 μ 放进 F_x 的计算里做闭环反馈，实际实现时你要选择：要么在状态方程里直接用 μ_hat 计算 F_x，让系统预测阶段就用当前估计的附着系数；要么把 F_x 当成外部扰动，只用可测制动压力换算。我的经验是前一种更稳健，因为它让滤波器自己内部形成了“μ 影响车轮减速 → 车速轮速差变化 → 观测残差反馈修正 μ”的闭合链路。

2.3 观测方程与传感器输入准备

车辆上常见的可测信号包括轮速、纵向加速度、横摆角速度。这个方案里观测向量取 z = [v_meas, ω_meas]。v 在实车上可以通过轮速换算基准加 GPS/惯导融合得到，仿真环境里就简单了，给真实信号叠加高斯白噪声即可。

观测方程是线性的：

h(x) = [v, ω]ᵀ

这么设计有个好处：EKF 里观测矩阵 H 是常数，可以少算一半 Jacobian，滤波器对角速度估计的依赖也更直接。

2.4 可观测性准备

在搭滤波之前建议先做个直观检查：系统状态是否可观测。v 和 ω 都能直接测，问题就集中在 μ 是否可以从测量序列中唯一反推。车轮制动时，纵向力改变 v 和 ω 的收敛速度，于是 μ 的影响通过 F_x 渗透到两个可测信号里，可观测性成立。而车轮接近纯滚动时 F_x 很小，v 和 ω 的运动主要由行驶阻力和输入力矩决定，μ 的影响微弱，可观测性退化。这一点后面调试时会专门讲到，也是我认为这个课题最需要提前理解的屏障。

3. EKF 与 UKF 的原理差异：雅可比矩阵、sigma 点与收敛表现

3.1 EKF：Jacobian 矩阵的那一步近似

EKF 的思路不复杂：对非线性系统做一阶线性化，然后用标准卡尔曼滤波的公式。预测方程：

x̂⁻(k+1) = f(x̂(k))

P⁻(k+1) = F·P(k)·Fᵀ + Q

F 是状态转移函数 f 对状态 x 的 Jacobian。更新方程：

K(k) = P⁻(k)·Hᵀ·(H·P⁻(k)·Hᵀ + R)⁻¹

x̂(k) = x̂⁻(k) + K(k)·(z(k) - h(x̂⁻(k)))

P(k) = (I - K(k)·H)·P⁻(k)

这个方案里 F 矩阵是 3×3，元素包括：

∂f_v/∂v、∂f_v/∂ω、∂f_v/∂μ、∂f_ω/∂v、∂f_ω/∂ω、∂f_ω/∂μ、∂f_μ/v 等。其中 μ 相关项来源于 F_x = μ·F_z 对 v、ω、μ 的导数，需要结合 Burckhardt 模型求链式导数。

EKF 的主要问题在强非线性场景。车速和轮速的动态与 μ 的关系在高滑移率区域并不是光滑线性，尤其路面类型突变时，一阶近似容易让协方差矩阵失真，估计值震荡甚至发散。但考虑计算成本，EKF 依然是嵌入式部署优先选择，而且在这个模型构架下，如果采样频率在 100Hz 以上，非线性影响通常可控。

3.2 UKF：用 sigma 点代替导数的确定性采样

UKF 的思路是“我不求导了，我直接让一群采样点穿过非线性函数”。对于 n 维状态，选取 2n+1 个 sigma 点：

χ₀ = x̄

χᵢ = x̄ + (√((n+λ)·P))ᵢ，i = 1...n

χᵢ = x̄ - (√((n+λ)·P))ᵢ，i = n+1...2n

权重为：

W⁰_m = λ/(n+λ)

W⁰_c = λ/(n+λ) + (1 - α² + β)

Wⁱ_m = Wⁱ_c = 1/(2·(n+λ))

参数 λ = α²(n + κ) - n，α 取小量如 1e-3，κ 通常取 0，β 对高斯分布取 2。

每个 sigma 点通过状态转移函数后，用加权统计重新得到预测均值和协方差，再对测量做同样的 sigma 点传播，计算交叉协方差得到卡尔曼增益。

UKF 不需要 Jacobian，这在系统里省掉很多手工推导，尤其当你后续换更复杂的轮胎模型（比如 Magic Formula）时，EKF 的 Jacobian 推导会变得非常痛苦，UKF 的优势一下就出来了。代价是计算量大约是 EKF 的 2-3 倍，因为每个时刻要传递 7 个 sigma 点（3 维状态）完成预测和更新两轮传播。

3.3 两相对比与工程取舍

我自己在同一套模型下跑过两者的对比，归纳进一张表：

我最终调试时用的是 UKF，因为它对初始化误差更宽容。初值 μ 设成 0.5，即使真实路面是压实雪，UKF 也能在 1 到 2 秒内收敛到 0.2 附近；EKF 在同样的初值误差下，第一阶段的协方差收缩过快，后续修正明显更慢。但如果你的目标是把代码烧进 MCU，EKF 其实也够用，要先用 UKF 把模型和调参逻辑理顺，再切换成 EKF 做量产代码。

4. Matlab Function 与 Simulink 模型搭建：从数据生成到滤波闭环

4.1 Simulink 顶层结构与信号流

我搭的 Simulink 模型整体分四个功能区。首先是激励输入模块，用于生成制动扭矩序列和路面 μ 序列，模拟一段先干沥青后冰雪的路面。其次是车辆动力学模块，用 Burckhardt 模型计算真实纵向力，积分得到真实车速和轮速，这部分相当于“真实系统”。然后是传感器模块，给 v 和 ω 叠加高斯白噪声，得到带噪声的量测。最后是估计器模块，用 UKF 或 EKF 从带噪声量测里估计 v、ω、μ，并和真实值做对比。

整个模型顶层信号流可以这样走：激励输入同时进入车辆动力学模块和估计器（因为制动力矩可以作为已知输入），车辆动力学模块输出真实 v、ω 给传感器模块，传感器模块输出带噪声测量给估计器，估计器输出估计状态给 Scope 做对比。

我用的是固定步长离散仿真，步长取 0.01 秒，对应 100Hz 更新频率。这里有个容易被忽视的点：真实车辆动力学模块的积分步长应该比估计器更新周期更小，否则你等于用一个粗离散系统去验证另一个粗离散系统。我习惯真实车辆动力学内部用变步长或 0.001 秒的固定步长，估计器每 10 步采集一次量测。

4.2 UKF 核心函数代码

这里给一个可直接跑的 UKF 主循环代码。状态 x = [v; omega; mu]，输入 u = [T_b; F_z]，量测 z = [v_meas; omega_meas]。为了简洁，状态转移函数 f 和观测函数 h 用嵌套函数写。代码路径基于标准的无迹变换，alpha、beta、kappa 按经验取值。

function [x_hat, P_est] = ukf_estimation(x_hat, P_est, z, u, dt, Q, R, params) % 参数 alpha = 1e-3; beta = 2; kappa = 0; n = numel(x_hat); lambda = alpha^2 * (n + kappa) - n; % sigma点生成 sqrtP = chol((n + lambda) * P_est, 'lower'); X = [x_hat, x_hat + sqrtP, x_hat - sqrtP]; % 权重 Wm = [lambda/(n+lambda), repmat(1/(2*(n+lambda)), 1, 2*n)]; Wc = Wm; Wc(1) = Wc(1) + (1 - alpha^2 + beta); % 预测：sigma点通过状态转移函数 n_points = size(X, 2); X_pred = zeros(n, n_points); for i = 1:n_points X_pred(:, i) = state_transition(X(:, i), u, dt, params); end x_pred = zeros(n, 1); for i = 1:n_points x_pred = x_pred + Wm(i) * X_pred(:, i); end P_pred = Q; for i = 1:n_points dX = X_pred(:, i) - x_pred; P_pred = P_pred + Wc(i) * (dX * dX'); end % 测量更新：sigma点通过观测函数 nz = numel(z); Z_pred = zeros(nz, n_points); for i = 1:n_points Z_pred(:, i) = observation_function(X_pred(:, i)); end z_pred = zeros(nz, 1); for i = 1:n_points z_pred = z_pred + Wm(i) * Z_pred(:, i); end Pzz = R; Pxz = zeros(n, nz); for i = 1:n_points 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 K = Pxz / Pzz; x_hat = x_pred + K * (z - z_pred); P_est = P_pred - K * Pzz * K'; end function x_next = state_transition(x, u, dt, params) v = x(1); omega = x(2); mu = x(3); T_b = u(1); F_z = u(2); R = params.R; I_w = params.I_w; m = params.m; % 纵向滑移率 if v > 0.1 slip = (v - omega * R) / v; slip = max(0, min(slip, 1)); else slip = 0; end F_x = mu * F_z; v_next = v + dt * (F_x / m); omega_next = omega + dt * (T_b - F_x * R) / I_w; mu_next = mu; % 随机游走 x_next = [v_next; omega_next; mu_next]; end function z_pred = observation_function(x) z_pred = [x(1); x(2)]; end

这个函数里有一个工程选择要说明：滑移率计算用了 max/min 限幅，防止 v 很小或轮速异常时滑移率出现负值或大于 1。在 v 接近 0 时直接置滑移率为 0，避免除零。

4.3 初始参数推荐与仿真步长

滤波器参数的初值我通常按下表来设：

如果你想验证 EKF，把上面函数换掉计算 Jacobian 那部分即可。EKF 里 F 矩阵三阶导数推导比较繁琐，但实现之后跑起来的曲线在强激励阶段和 UKF 几乎重合，差异主要体现在路面切换瞬间和弱激励阶段。

4.4 Simulink 里 MATLAB Function 模块的具体搭建步骤

在实际 Simulink 模型里，估计器我是用 MATLAB Function 模块实现的。创建步骤如下：

1. 从模块库拖入 MATLAB Function 模块，打开代码编辑器。
2. 将上面 ukf_estimation 函数主体完整粘贴，把参数结构体 params 改为单独输入端口（R、I_w、m），避免结构体在 Simulink 代码生成时出问题。
3. 模块端口设计为四输入两输出：输入为量测 z、制动扭矩 T_b、法向载荷 F_z、采样周期 dt，输出估计 x_hat 和协方差 P_est。
4. 在 MATLAB Function 模块内部用 persistent 变量保存 x_hat、P_est，保证跨采样步的状态保持。
5. 模块内部每次调用时先取 persistent 值，调用 ukf_estimation，再把新 x_hat、P_est 存回 persistent。

persistent 变量的使用有个小坑：Simulink 默认会在第一个仿真步把 persistent 初始化成空矩阵，所以代码里必须判断 isempty 做初始化。我习惯把初值定义在 persistent 变量之外的一个常量结构体里，比如 Editor PreLoadFcn 回调。

模型跑通后，建议先不开噪声，让滤波器去跟踪理想信号，确认状态转移和更新逻辑没有 bug。跟踪误差在滑移率快速变化阶段会有一个小凸起，那是滤波器的固有响应延迟，不是程序问题。再叠加噪声后，观察估计曲线是否有明显滞后或发散。

5. 仿真调试中的观测性问题、噪声标定与收敛性陷阱

5.1 车轮接近纯滚动时的可观测性退化：这个坑最隐蔽

一旦进入小滑移率区间（比如匀速巡航），车轮和路面之间几乎没有纵向力激励，μ 对观测信号的影响趋近于零。此时协方差矩阵中 μ 对应的方差会慢慢变大，滤波器的修正增益趋近于一个激进的状态，一旦量测噪声带来一个尖峰，估计值可能瞬间跳变。

我在仿真里遇到过这样的情况：车辆先以 20m/s 匀速跑了 3 秒，然后才开始制动。UKF 已经收敛到 μ=0.8 附近，匀速期间虽然没有新的信息进来，但 μ 估计值居然慢慢漂到了 0.95 以上。原因就是随机游走的 μ 状态在无激励时只能靠过程噪声维持方差，更新环节几乎不修正它，漂移源于数值误差和协方差膨胀。

对策有三种。第一是给 μ 的过程噪声调小，让它不至于在无激励时乱跳，但同时也会让真正路面切换时收敛变慢。第二种是我推荐的：引入激励检测逻辑，只有当滑移率绝对值的统计量超过某一阈值（比如连续 20ms 大于 0.05）才允许 μ 维度的增益显著起作用。第三种是在状态方程里给 μ 加一个回拉项，比如 μ_next = μ + k·(μ_ref - μ)，μ_ref 取一个保守的中值，但这会把估计结果向先验拉偏，跟踪路况变化的能力会下降。实际项目中激励检测用得最多。

5.2 Q/R 标定节奏：别一上来就调六个参数

新手拿到滤波问题最容易干的事是六个参数来回试，最后越调越晕。我总结的节奏是三步走。

第一步，固定 R，把 Q(1,1) 和 Q(2,2) 调到一个能让 v 和 ω 的估计曲线与量测曲线足够接近、但又不至于完全跟噪声走的水平。标准是估计曲线的“平滑度”和真实曲线的“趋势保真度”都要合格。车速和轮速是直接可测的，这个阶段状态 x 的估计精度不依赖 μ。

第二步，保持 Q(1,1)、Q(2,2) 不变，单独调 Q(3,3)。把 μ 的初值故意设偏，比如真实 0.8 设成 0.3，观察收敛时间和超调量。收敛时间在 1 秒以内算快，超调太大会导致 v 估计短暂跟着跑偏，这时候就该降低 Q(3,3)。如果收敛太慢但最终能到，就增大 Q(3,3)。

第三步，把量测噪声 R 调准。R 应该对应传感器实测噪声方差，不要拿仿真里的噪声方差直接填，因为实际传感器还有量化噪声和安装误差。最稳妥的方法是采集一段车辆静止时的轮速信号，计算方差作为 R 的参考值。

调参指标我建议配合两个评价量：估计值和真值的均方根误差（RMSE），以及 μ 估计从初值到真值附近的上升时间。RMSE 反映稳态精度，上升时间反映动态响应。同时也要观察协方差矩阵对角元的合理性，P(3,3) 如果长时间过大，说明激励不够或 Q 设置偏大。

5.3 量测离群值与延时处理

实车信号不会像仿真这么干净。轮速传感器偶尔会丢脉冲，CAN 信号偶发跳变，这些离群值对卡尔曼滤波的伤害比高斯噪声大得多，因为它会让增益矩阵瞬间放大，把状态估计拉偏后还要好几步才能缓过来。

比较实用的保护是在更新之前做一个残差门限判断。计算新息向量 v_new = z - h(x_pred)，然后检查其范数是否超过根据 Pzz 设定的阈值，比如 3 倍标准差。超了就直接跳过本次量测更新，让预测值往下走。代价是离群值持续发生时滤波器会退化成纯预测，但这在真实底盘控制器里比干脆发散掉强得多。

对于信号延时的问题，可以在模型中加一个补偿：如果量测延迟已知为 τ，那么比较量测和预测时，把预测状态也回退 τ 产生对应的预测量测。但这个方案会增加复杂度，项目初期可以先不管延迟，仿真步长 0.01s 时量测延迟一两步的影响很小。

5.4 路面切换瞬间的表现与激励需求

路面从干沥青切到冰雪时，真实 μ 跳变，但滤波器中的 μ 是慢变量，不可能瞬间跟上。这里有个值得注意的现象：路面切换后，如果车辆仍在稳定制动，滑移率还在峰值左侧，那么 μ 只能以一个相对缓慢的速度爬向新值。想让收敛快一点，必须让车轮产生更大的滑移激励，也就是瞬时降低制动压力或转向扰动，把轮胎推进非线性区。这部分依赖上层控制策略的配合，不是滤波器本身能完全决定的。

在仿真模型中，我习惯在路面切换点附近给制动扭矩增加一个 200ms 的方波脉冲，模拟驾驶员急刹，以便观察滤波器在最恶劣条件下的跟踪能力。用这样一组仿真数据来评估 UKF 和 EKF 的差异最客观，否则很可能会出现“看起来差不多”的误导结论。

我在实际项目里通常保留两套估计器：算法验证阶段全部用 UKF，占个“模型改动不用重推导数”的便宜；到了要生成 C 代码嵌入控制器的阶段，再换成 EKF，把 UKF 阶段调好的噪声参数直接迁移过去，重新跑一轮仿真对比。如果你也准备走这个路线，我建议在 Simulink 里把估计器封装成独立的子系统，输入输出接口固定，两套算法只换内部实现，模型上层结构完全不用动。这样切换成本最低，也方便回归对比。