路面附着系数估计实战：基于EKF与UKF的Simulink实现与对比-北京尧图网络

做车辆动力学控制方向的朋友几乎迟早都会撞上同一个问题当前路面到底能提供多大的附着系数。ABS要靠它判断车轮会不会抱死ESP要靠它决定要不要介入AEB在低附着路面上能不能在目标距离内刹停也跟它直接挂钩。这个参数在整车算法里就是个关键前提但麻烦在于它没法直接装传感器测出来只能靠估计。这篇文章要聊的就是我在Matlab/Simulink里完成的一个路面附着系数估计项目分别用EKF和UKF两种滤波方案做了实现和对比从车辆动力学建模、轮胎模型、滤波器设计到仿真调试整条链路完整走了一遍过程中踩过的坑也一并分享出来。适合正在做车辆状态估计课题的研究生、刚接触ADAS算法的工程师参考。先说明一下项目背景。整车控制器里很多功能模块都需要μ值但μ本身受路面材料、轮胎磨损、车速、滑移率等多种因素影响是典型的时变参数。直接测量需要昂贵的光学传感器或专用测试轮胎量产车上不可能装。工程上的常规做法是基于可测信号比如轮速、纵向加速度、转向角通过车辆动力学的约束关系把μ“反演”出来。这个反演过程本质上就是一个非线性状态估计问题所以核心工具自然落到了非线性卡尔曼滤波上。我用Simulink搭了一套完整的仿真验证环境包含车辆单轮动力学模型、魔术公式轮胎模型、EKF/UKF估计器模块和传感器噪声注入模块通过切换不同附着系数的路面工况对比了两种滤波器的估计精度、收敛速度和鲁棒性。整体下来UKF在强非线性工况下的优势很明显但EKF的计算开销确实小不少两者如何取舍后面我会详细展开。1. 算法选型EKF和UKF为什么不是普通卡尔曼滤波1.1 从KF到EKF线性化到底引入了什么问题经典卡尔曼滤波KF要求系统是线性的也就是状态方程和观测方程都可以写成矩阵乘法的形式。但车辆动力学系统天然是非线性的轮胎力跟滑移率之间是那种带饱和特性的曲线关系稍微算重一点就是强非线性。这时候第一个想到的扩展办法就是扩展卡尔曼滤波。EKF的思路很直接在每一个时刻的估计值附近把非线性函数做一阶泰勒展开用雅可比矩阵替代原来的状态转移矩阵和观测矩阵然后套用标准KF的预测-更新流程。工程上用起来很顺手因为计算量跟KF基本一个量级很多嵌入式平台都能跑。但EKF有两个明显的痛点。第一雅可比矩阵求起来麻烦尤其是轮胎模型这种嵌套了很多三角函数的表达式手推偏导容易错用符号工具箱生成又费时间。第二一阶线性化会丢高阶项在强非线性区域误差会被放大甚至导致滤波器发散。实际调试中我遇到过几次估计值突然跳变的状况排查下来就是线性化误差惹的祸。1.2 UKF用采样代替线性化无迹卡尔曼滤波UKF走了另一条路不强行线性化而是用无迹变换来传播状态分布。它的核心思路是选取一组确定的Sigma点让这些点的均值和协方差等于当前状态估计的均值和协方差然后把这组点一个一个送入非线性函数通过输出的统计量来近似真实的后验分布。这样做的好处很明显不用求雅可比矩阵实现代码里面少掉一大块容易出错的部分精度上能保留到二阶矩比EKF的一阶展开高。代价就是多了一次Sigma点的传播计算状态维度越高计算量增长越快。在我这个三状态模型里多出来的开销其实可以接受。1.3 实际项目里怎么选选EKF还是UKF不能只看精度得结合项目约束来权衡。如果控制器是量产ECU算力受限、代码要过MISRA规范EKF相对好落地雅可比矩阵虽然推导麻烦但可以提前求好解析表达式固化下来。如果是在做预研、算法验证或者课题研究UKF是非常好的选择代码结构清晰调试周期短出问题的概率低。我做这个项目时两边都实现了因为想对比一下在路面附着系数估计这个具体场景下两者的差距到底有多大。结论先说在干沥青和湿沥青这类常规路面上两者差异不大但在冰雪路面这种极端低附着、轮胎力进入饱和区的工况下UKF的估计稳定性明显好于EKF。后面章节我会给出具体的仿真对比。2. 车辆动力学建模估计器的“心脏”2.1 单轮车辆模型估计器本身只是一个数学算法它要知道车辆的运动规律才能把状态“猜”出来。所以我先建立了简化的车辆纵向动力学模型。用单轮模型的原因是这个项目聚焦纵向工况不考虑转向和侧向力把问题约束在一条直线上模型简单且能覆盖绝大多数ABS/AEB场景。单轮模型的核心是两个方程。第一个是整车纵向运动方程m * dv/dt -F_xm是整车质量的四分之一v是车身纵向速度。第二个是车轮旋转运动方程J * dω/dt T_drive - F_x * R驱动工况 J * dω/dt -T_brake F_x * R制动工况其中J是车轮转动惯量ω是车轮角速度R是轮胎滚动半径T_drive和T_brake是作用在车轮上的驱动和制动力矩。这两个方程建立起车身速度、轮速和轮胎力之间的动态关系。轮胎力F_x是连接两者的关键变量而F_x又由路面附着系数μ决定。这就构成了状态估计的基本框架给滤波器一个多自由度系统输入端是可以操控的力矩可测量的输出是轮速和车身速度需要估出来的是隐含的μ。2.2 魔术公式轮胎模型轮胎力跟滑移率之间的关系行业内最常用的是Pacejka提出的魔术公式。它用一组三角函数拟合轮胎实验数据表达式如下F_x D * sin(C * arctan(B * κ - E * (B * κ - arctan(B * κ))))其中κ是滑移率B是刚度因子C是形状因子D是峰值因子E是曲率因子。D值在峰值附着系数附近直接决定轮胎力的峰值工程上近似取D μ * F_z所以路面附着系数μ的变化会直接体现为轮胎力曲线的整体缩放。这里我补充一下滑移率的定义。制动工况下κ (v - ωR)/v驱动工况下κ (ωR - v)/(ωR)取值范围大致在0到1之间。可以这样理解滑移率表示车轮“打滑”的程度纯滚动时κ接近0车轮完全抱死时κ 1。轮胎力跟κ的关系是一条先上升后饱和的曲线在κ 0.1到0.2附近轮胎力达到峰值之后进入饱和区。魔术公式的好处是参数少、精度高而且对μ的变化非常敏感这正是估计器需要的特性。模型参数我参考了典型的乘用车轮胎数据参数数值说明B10刚度因子C1.3形状因子Dμ * F_z峰值因子E0.97曲率因子F_z3460 N单轮垂直载荷R0.3 m轮胎滚动半径J0.6 kg·m²车轮转动惯量m353 kg四分之一整车质量2.3 状态方程和观测方程有了模型接下来就是推导状态方程。我把状态向量选为x [v, ω, μ]^T把μ当作一个慢时变参数扩展进状态里用随机游走模型描述它的变化μ_dot w_μw_μ是过程噪声。连续时间状态方程写出来就是这样dv/dt -μ * F_z * sin(C * arctan(B * κ - E * (B * κ - arctan(B * κ)))) / mdω/dt (T_brake μ * F_z * R * sin(C * arctan(B * κ - E * (B * κ - arctan(B * κ)))) / Jdμ/dt 0注意这里我把F_x用魔术公式展开κ又是v和ω的函数所以状态方程整体是v、ω、μ的强非线性耦合方程组。观测方程方面轮速传感器可以直接测ω如果车辆装有纵向加速度传感器还能测到a_x -F_x/m。所以我设观测向量y [ω, a_x]^T。写成观测方程y1 ω y2 -μ * F_z * sin(C * arctan(B * κ - E * (B * κ - arctan(B * κ)))) / m这里有个非常关键的工程细节如果没有加速度传感器单靠轮速信号μ的观测性会很差。仿真中你会发现附着系数估计值迟迟不收敛或者收敛得很慢。加入加速度信号后相当于给估计器增加了一个直接观测轮胎力的通道μ的可观测性大大改善。实际项目中如果只有轮速信号就需要额外的激励比如制动或加速过程来激发可观性。3. Simulink实现从模型到仿真环境3.1 整体架构设计我把整个Simulink模型拆成了几个层次分明的模块方便独立调试和替换算法。顶层主要是信号源、被控对象模型、估计器和结果显示四个部分。信号源产生制动力矩输入这里我用了阶跃信号和斜坡信号两种模拟驾驶员突然踩刹车和缓慢制动两种场景。被控对象模型严格按上一节的状态方程搭建输入力矩输出真实的车身速度、轮速和轮胎力。估计器模块接收带噪声的传感器信号和力矩信号实时输出μ的估计值。结果展示用Scope和To Workspace模块把估计值和真实值放到同一张图里对比。这个架构最大的好处是在仿真里可以拿到“真实值”。真实的车速、真实的μ在实车上永远测不准但仿真里模型自己知道这样就能非常直观地评估滤波器的估计误差。这个优势一定要利用好后面的参数调优全靠它。3.2 EKF/UKF模块搭建细节估计器我用了Simulink的MATLAB Function模块来实现。选择MATLAB Function而不是手写S-Function主要是考虑到代码可读性和开发效率。MATLAB Function可以在Simulink里直接调试也能用MATLAB的调试器设置断点对排查算法逻辑问题太方便了。MATLAB Function模块的输入端口我设计了四个轮速测量值、加速度测量值、制动力矩和采样时间。输出端口两个估计的μ值和估计的车身速度。模块内部需要保存上一时刻的状态和协方差矩阵所以在函数内部用persistent变量来存储初始化时赋单位矩阵或零矩阵。这里有个细节提醒一下persistent变量在Simulink模型运行前不会被自动重置。如果你改了算法参数重新跑仿真最好先执行一次clear all或者在模型初始化回调里加复位语句不然上一轮的残留值会污染新一轮仿真的初值。这个坑我踩过一次调试时发现第二次运行的结果跟第一次不一样排查了半天才知道是persistent变量没清。3.3 EKF核心代码实现EKF实现的核心是按predict和update两步走。predict步骤里我用四阶龙格库塔积分对状态方程做离散化把连续时间模型转为离散递推形式。这一步比简单的欧拉法精度高不少但计算量也大一些在Simulink仿真里完全没问题。EKF算法关键流程我写成了一段MATLAB代码供参考function [mu_est, v_est] ekf_estimator(w_meas, ax_meas, T_brake, dt) % 状态: x [v, omega, mu] persistent x P Q R if isempty(x) x [20; 60; 0.8]; % 初始速度20 m/s, 轮速60 rad/s, 初始mu0.8 P eye(3) * 0.1; Q diag([0.5, 2, 0.01]); % 过程噪声mu的噪声要设小 R diag([0.1, 0.5]); % 测量噪声轮速和加速度传感器 end % 1. 预测步使用四阶龙格库塔 x_pred rk4_state(x, T_brake, dt); % 计算状态转移雅可比矩阵 F F calc_F_jacobian(x, T_brake); P_pred F * P * F Q; % 2. 更新步 % 计算观测雅可比矩阵 H H calc_H_jacobian(x_pred); y_meas [w_meas; ax_meas]; y_pred [x_pred(2); calc_ax(x_pred)]; S H * P_pred * H R; K P_pred * H / S; x_upd x_pred K * (y_meas - y_pred); P_upd (eye(3) - K * H) * P_pred; % 保存状态并输出 x x_upd; P P_upd; mu_est x(3); v_est x(1); end这段代码的关键在雅可比矩阵计算。calc_F_jacobian需要对状态方程里三个状态分别求偏导calc_H_jacobian需要对观测方程求偏导。因为魔术公式里有嵌套的arctan和sin偏导表达式很冗长。我建议用MATLAB的Symbolic Math Toolbox先推导并验证再手动转成数值函数能大幅减少出错概率。3.4 UKF核心代码实现UKF实现起来反而比EKF更简洁因为不需要雅可比矩阵。核心是生成Sigma点、通过非线性函数传播、计算权重合并结果三个步骤。Sigma点个数是2n1个n是状态维度我这个项目n3所以每次递推要传播7个点。UKF算法关键流程代码function [mu_est, v_est] ukf_estimator(w_meas, ax_meas, T_brake, dt) % 状态: x [v, omega, mu] persistent x P Q R alpha beta kappa if isempty(x) x [20; 60; 0.8]; P eye(3) * 0.1; Q diag([0.5, 2, 0.01]); R diag([0.1, 0.5]); alpha 1e-3; % Sigma点分布参数 beta 2; % 高斯分布最优值 kappa 0; end n numel(x); lambda alpha^2 * (n kappa) - n; % 计算Sigma点 [U, S, ~] svd(P); sqrtP U * diag(sqrt(diag(S))) * U; X zeros(n, 2*n1); X(:, 1) x; for i 1:n X(:, i1) x sqrt(n lambda) * sqrtP(:, i); X(:, ni1) x - sqrt(n lambda) * sqrtP(:, i); end % Sigma点传播 Y zeros(n, 2*n1); for i 1:2*n1 Y(:, i) rk4_state(X(:, i), T_brake, dt); end % 权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); Wm(2:end) 1 / (2 * (n lambda)); Wc(2:end) 1 / (2 * (n lambda)); % 合并预测 x_pred sum(Y * Wm, 2); P_pred Q; for i 1:2*n1 diff Y(:, i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end P_pred (P_pred P_pred) / 2; % 保证对称性 % 更新步 Z zeros(2, 2*n1); for i 1:2*n1 Z(:, i) [Y(2, i); calc_ax(Y(:, i))]; end z_pred sum(Z * Wm, 2); Pzz R; for i 1:2*n1 dz Z(:, i) - z_pred; Pzz Pzz Wc(i) * (dz * dz); end Pxz zeros(n, 2); for i 1:2*n1 dx Y(:, i) - x_pred; dz Z(:, i) - z_pred; Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; y_meas [w_meas; ax_meas]; x_upd x_pred K * (y_meas - z_pred); P_upd P_pred - K * Pzz * K; P_upd (P_upd P_upd) / 2; x x_upd; P P_upd; mu_est x(3); v_est x(1); endUKF代码里有一个关键点计算协方差矩阵平方根时不能直接用Cholesky分解因为数值误差可能导致矩阵非正定。我用SVD分解替代鲁棒性更好。另外每次更新后对协方差矩阵做对称化处理避免数值误差累积导致矩阵不对称。4. 实验工况设计与参数调优4.1 典型路面工况设置仿真验证阶段我设计了三种典型工况来考察估计器的性能。工况一是干沥青路面μ真实值为0.85这是最常见的高附着场景。工况二是湿沥青路面μ真实值为0.5。工况三是冰雪路面μ真实值为0.2这是对估计器挑战最大的场景。每种工况下我都模拟了两种制动力矩输入。第一种是阶跃制动驾驶员在t2秒突然踩下制动踏板力矩从0瞬间跳到800 N·m模拟紧急制动。第二种是斜坡制动力矩从0在5秒内线性增加到800 N·m模拟渐进制动。为什么两种激励都要测因为路面附着系数估计的收敛速度和可观测性跟激励强度密切相关。阶跃激励能让轮胎迅速进入大滑移区域μ的信息更容易被激发出来斜坡激励更接近日常驾驶估计器需要在相对平缓的激励下也能维持稳定。两套测试合在一起才能完整评估估计器在不同激烈程度工况下的表现。4.2 噪声参数与滤波器参数标定滤波器参数标定是整个过程最花时间的一环。Q矩阵反映对模型不确定度的信任程度R矩阵反映对传感器测量值的信任程度。Q设得太大估计值会剧烈跳动跟随噪声Q设得太小滤波器反应迟钝真实状态变化时跟踪不上。我最终的参数设置如下参数数值设置依据Q(1,1)0.5对应车速过程噪声考虑模型简化的误差Q(2,2)2.0对应轮速过程噪声制动力矩波动影响Q(3,3)0.01对应μ的过程噪声μ变化慢噪声要小R(1,1)0.1轮速传感器噪声功率R(2,2)0.5加速度传感器噪声功率初始P0.1 * I对初值有一定不确定度μ对应的Q(3,3)是整个调参的枢纽。设太大μ估计值会跟随机噪声乱跑出现无物理意义的剧烈波动设太小路面变化时μ估计收敛极慢比如从干沥青开到冰雪路面可能要好几秒才能跟上这在安全场景中是不可接受的。我反复试下来取0.01到0.05之间比较合适。4.3 仿真结果EKF与UKF对比直接说结论。在干沥青工况下EKF和UKF的表现差距不大两者都能在大约1秒内收敛到真实值0.85附近稳态误差都在0.03以内。原因在于干沥青路面的轮胎力-滑移率曲线线性段较长非线性程度不高EKF的一阶线性化损失不明显。差距在冰雪路面工况下彻底拉开。冰雪路面μ0.2轮胎力在很小的滑移率下就进入饱和区非线性非常强。EKF的估计结果出现了明显的振荡最大误差达到0.08而且收敛后有持续的波动。UKF的估计曲线则平滑得多收敛后误差控制在0.03以内振荡幅度只有EKF的一半左右。这个结论其实跟学术文献里的结论一致轮胎力饱和区域是典型的强非线性区域EKF丢失的高阶信息恰恰是估计的关键。所以在做低附着路面估计时我强烈建议优先上UKF。5. 常见问题与排查技巧实录5.1 估计结果发散这是碰到最多的一个问题。现象是仿真跑着跑着μ估计值突然跳到离谱的数值比如变成负数或者超过2然后一骑绝尘不回头。排查下来常见原因有三个Q矩阵设置过大、协方差矩阵失去正定性、以及激励不足导致的可观测性缺失。激励不足的情况最隐蔽。如果车速和轮速信号基本恒定轮胎没有进入明显滑移状态μ的信息从测量中根本提取不出来估计器可能收敛到任意值。这时候解决方法是在仿真中加入足够的制动力矩激励让轮胎经历从弹性区到饱和区的完整过程。实车标定时同理纯匀速巡航时不要指望μ估计值有多准。5.2 收敛速度太慢μ的收敛慢多半是Q(3,3)设得太小了。我曾经为了抑制估计值的波动把Q(3,3)压到0.001结果在仿真中从干沥青切换到湿沥青路面时μ估计值花了将近5秒才跟上真实值这个响应速度在安全控制场景中完全不可用。后来调回0.01收敛时间压缩到了1.5秒左右。另外滤波初值也会影响收敛速度。设置一个合理的初始值能明显缩短收敛时间我一般按当前环境先猜一个μ比如晴天默认0.9雨天默认0.5而不是盲目用0或者1的默认值。5.3 数值稳定性问题EKF在运行时偶尔会遇到协方差矩阵非正定的情况导致Kalman增益计算出现NaN。根因通常是数值误差累积。我处理的办法第一个是每次更新后强制对称化用P (PP)/2第二个是在计算卡尔曼增益时改用稳健的求解方式第三个是给P的对角元素加一个很小的下界避免出现负方差。我在UKF的Sigma点生成中就遇到过这种情况直接用Cholesky分解报错改成SVD分解后问题消失。这类数值问题在纯数学推导中不会出现但到代码实现里天天见。建议你把数值稳定性检查写进代码的公共工具函数里可以少掉很多头发。5.4 问题速查表问题现象可能原因排查方向μ估计值发散到异常范围Q矩阵过大 / 协方差矩阵失去正定性调小Q检查P矩阵数值稳定性μ估计收敛过慢Q(3,3)过小 / 初值偏差过大调大Q(3,3)设置合理初值EKF和UKF差距不大工况在线性区 / 激励不足加大制动力矩激励UKF估计曲线有毛刺过程噪声参数偏大 / 采样时间过长调小Q减小仿真步长二次仿真结果与首次不同persistent变量未清零运行前执行clear all6. 从仿真到实车的几点心得仿真环境跑通只是第一步真要往实车上迁移有几个点需要提前想清楚。实车传感器信号不像仿真里给的那么干净轮速传感器在低速时分辨率变差加速度传感器有安装偏置这些都需要在算法外面加信号处理。我的建议是设计一个前置的数据有效性判断模块在车速过低或者传感器自检异常时暂停估计器更新冻结上次输出。关于算力EKF在2GHz的工控机上单步运行耗时约0.15毫秒UKF大约0.4毫秒。如果控制器主频不高或者有其他算法共享CPU这个差异还是值得考虑的。我的个人经验是如果目标是量产项目先评估一下系统非线性强度如果大部分工况都在轮胎曲线的线性段附近EKF完全够用如果追求极限工况的性能或者不想推导雅可比矩阵直接上UKF。最后再分享一个细节仿真中调参找到一个合适的Q和R之后实车调试时往往需要把R适当调大。因为仿真里的传感器噪声是仿真器给的参数完全可知实车的噪声特性比较复杂存在未建模的干扰把R调大一些相当于承认“测量没那么可信”反而有助于算法稳定。我在台架测试中就吃过这个亏仿真里调好的参数上车直接振荡把R调大一倍之后就好了。这个项目做完之后我最大的感受是路面附着系数估计的难点真的不在滤波算法本身而在于你愿不愿意把车辆动力学模型的机理搞清楚。模型建准了EKF和UKF都差不到哪去模型建歪了再高级的滤波算法也救不回来。希望这篇文章能帮你在Simulink里把整个链路跑通少走一些我走过的弯路。

路面附着系数估计实战：基于EKF与UKF的Simulink实现与对比

相关推荐

阿兹海默症MRI辅助诊断系统：从CNN到Grad-CAM的深度学习实践

VisualSVN Server+TortoiseSVN从零搭建SVN版本控制系统完整指南

LPC1114例程详解：从GPIO到外设驱动开发实战

更多资讯

Turso 异步 I/O 模型深入解析：协作式让出、显式状态机与 CompletionGroup 实践

界面开发1.0实战：从零搭建完整前端项目的经验与坑位复盘

企业微信API与RPA实现CRM自动化对接实战

spaCy 如何用 spacy benchmark speed 测量流水线吞吐量？

树莓派Pico实战：电位器控制LED的ADC与PWM应用

OpenClaw+优云智算+Coding Plan：实现从灵感到发布的全流程自动化

今日更新

SpacetimeDB C++ 模块库实战：用 C++20 在数据库内构建 WebAssembly 模块

别踩雷！不是所有 AI 写作工具都靠谱，2026 高校认可工具精选

SSM+MySQL幼儿园管理系统源码解析：从部署到优化

本周精选

超人会飞不算本事：系统稳定依赖清晰规则与边界设计

超人VS蜘蛛侠：拆解超级IP的影响力与传播方法论

基于CNN的调制信号识别：MATLAB实现时频图分类实战

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？