CN103407451B - 一种道路纵向附着系数估计方法 - Google Patents

Abstract

本发明公开了一种道路纵向附着系数估计方法。在平坦的高速公路路面上，针对前轮转向四轮汽车，在整车纵向动力学模型和简化魔术公式轮胎模型的基础上，利用带遗忘因子的递归最小二乘（Recursive least squares,RLS）方法实时初步估计出纵向道路附着系数，进一步将所估计出的纵向道路附着系数与轮胎模型参数作为扩充状态，利用扩展卡尔曼滤波（Extended Kalman Filter,EKF）算法，滤除信号噪声，并实现轮胎模型系数的自适应调整，最终实时获取准确、鲁棒的道路纵向附着系数估计，该方法能够同时适应平坦高速公路上的高滑移率和低滑移率工况。

Description

一种道路纵向附着系数估计方法

技术领域

本发明涉及一种道路纵向附着系数估计方法，其目的在于通过整车纵向动力学建模和轮胎模型，利用递归最小二乘和扩展卡尔曼滤波方法，在平坦的高速公路路面环境下，实现在高滑移率和低滑移率工况下对纵向道路附着系数的准确、可靠估计，汽车纵向主动安全***可以根据所估计出的纵向道路附着系数调节控制策略，提高车辆安全，具有精度高、成本低、实时性和适应性好等显著优点，属于汽车主动安全测量及控制领域。

背景技术

随着社会经济的发展，道路交通安全问题日益突出，并已成为全球性难题。全世界每年因交通事故都会造成大量的人员伤亡和财产损失，世界各国都在努力降低交通事故的发生。近年来，汽车主动安全技术得到了迅速的发展。汽车主动安全技术能够防患于未然，主动避免事故的发生，已成为现代汽车最主要的发展方向之一。目前常见的主动安全技术主要包括防抱死制动***(ABS)，车辆电子稳定程序(ESP)，牵引力控制***(TCS)，电控驱动防滑***(ASR)，四轮转向稳定控制***(4WS)等。这些***效果的优劣很大程度上取决于能否“道路自适应”，即如果能够实时估计出道路附着系数，***就可以根据当前路况调节控制策略，提高车辆安全。纵向道路附着系数是影响汽车制动性能的一个关键因素，是碰撞预警／主动避撞(CW／CA)、自适应巡航(ACC)、防抱死制动***(ABS)等汽车纵向安全辅助***中必不可少的一个重要参数，实时估计纵向道路附着系数是设计符合实际要求的纵向安全算法的重要技术前提，其精度直接关系汽车的行驶安全性与稳定性，即纵向主动安全控制***能否有效工作在很大程度上依赖于道路附着系数能否被实时、准确的测量或估计。

目前，在汽车主动安全领域，道路附着系数主要分为直接测量和间接估计两类，直接测量方法是利用光、声、图像、雷达等传感器直接检测路面，测量一些对路面附着系数影响较大的因素，并根据以往经验预测当前道路附着系数的大小，但这些方法都需要额外加装传感器，且传感器成本都较高，难以实现大规模的商业应用；其次需要进行大量的测试训练，识别精度很大程度上依赖于经验，难以准确估算没有测试和训练过的路况的附着系数。间接估计方法是通过对汽车的运行过程进行运动学或动力学建模，结合轮胎模型，将有关低成本的车载传感器(如轮速传感器、陀螺仪、加速度计以及GPS等)信息作为观测信息，进而利用适当的滤波估计算法实现对道路附着系数的估计。已有的道路附着系数间接估计方法包括基于车辆侧向动力学和基于纵向动力学的研究两种，其中，前者主要用于汽车转向运动，而对于汽车纵向安全辅助***，较为关注的是基于车辆纵向动力学的附着系数估算。不同附着系数的路面上具有不同的滑移率和归一化牵引力关系，这是利用车辆纵向动力学估算附着系数的基础。基于这种关系，就可以通过估计滑移斜率，利用分类算法来估计附着系数，这是当前国内外学者所使用最多的方法，该方法认为小滑移率范围内的归一化牵引力与滑移率成正比，其比例关系称为滑移斜率，而道路附着系数与滑移斜率成正比，因此，通过估计滑移斜率，就能估计出轮胎路面的附着系数。但这种方法只能应用于低滑移率工况，而在高滑移率工况下，归一化牵引力与滑移斜率成正比的假设不再成立，因此无法利用该方法进行道路附着系数估计；同时在同一路面上滑移斜率随其他因素的影响也存在一定差异，给附着系数的识别带来了一定的困难。另一种常用的方法是采用不同形式的非线性公式实时根据测量数据拟合归一化牵引力与滑移率曲线，从而根据计算所拟合曲线的峰值估计附着系数，但这种方法需要在较大轮胎滑移率条件下才能使估算结果更加接近于真实值。同时，该方法采用非线性函数进行拟合，计算量大，实时性较差，且不论采用哪一种非线性函数，都无法在所有工况下都吻合轮胎特性，误差较大。可见，当前在附着系数估计上所存在的主要问题在于难以同时在高滑移率和低滑移率工况下都取得较好的估计效果，限制了其在汽车纵向主动安全***上的应用。

发明内容

为实现在高滑移率和低滑移率工况下都能对道路附着系数进行准确、可靠估计，本发明提出了一种道路纵向附着系数估计方法。本发明针对现有方法难以同时在高滑移率和低滑移率工况下都取得较好的估计效果，面向平坦的高速公路上所行驶的前轮转向四轮汽车，采用更符合实际车辆运行过程的简化魔术公式轮胎模型，结合纵向整车动力学模型利用递归最小二乘法进行道路附着系数估计，同时充分利用低成本车载传感器信息作为***的观测量，通过扩展卡尔曼滤波方法实现轮胎模型系数的自适应调整，并去除纵向附着系数估计结果所含噪声，输出最终的道路纵向附着系数估计结果，具有精度高、成本低、实时性和适应性好等显著优点。

一种道路纵向附着系数估计方法，其特征在于：在平坦的高速公路路面上，针对前轮转向四轮汽车，基于整车纵向动力学模型和简化魔术公式轮胎模型，利用带遗忘因子的递归最小二乘(Recursive least squares，RLS)方法实时初步估计出道路纵向附着系数，进一步将所估计出的道路纵向附着系数与轮胎模型参数作为扩充状态，利用扩展卡尔曼滤波(Extended Kalman Filter，EKF)算法，滤除信号噪声，并实现轮胎模型系数的自适应调整，最终实时获取准确、鲁棒的道路纵向附着系数估计，能够同时适应平坦高速公路上的高滑移率和低滑移率工况；

具体步骤包括：

1)计算纵向滑移率

用s_j(j＝f，r)表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_f和后轮轴纵向滑移率s_r，下角标j取f或r，f或r分别表示前或后轮轴，s_j计算方法为：

式(1)中，R_e表示轮胎半径；v_tf和v_tr分别表示前、后轮轴上沿轮胎方向的速度，本发明中都以车辆绝对速度v代替，通过GPS测量获取；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_j(j＝f，r)且

${\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{fR}}+{\mathrm{\ω}}_{\mathrm{fL}})$

                           (2)

${\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{rR}}+{\mathrm{\ω}}_{\mathrm{rL}})$

式(2)中，ω_fL、ω_fR、ω_rL和ω_rR分别表示左前轮、右前轮、左后轮和右后轮的旋转角速度，通过利用四个轮速传感器测量获得；

2)计算整车纵向力

用表示归一化牵引力，即又可分为前轮轴归一化牵引力和后轮轴归一化牵引力下角标j取f或r，f或r分别表示前或后轮轴：

式(3)中，F_xf和F_xr分别表示前、后轮轴上纵向力，F_xf和F_xr可统一记为F_xj(j＝f，r)，F_zf和F_zr分别表示分配到前或后轮轴上的垂向载荷，可统一记为F_zj(j＝f，r)且可按下式计算：

$F_{\mathrm{zf}}=\frac{\mathrm{mgb}}{(a+b)},F_{\mathrm{zr}}=\frac{\mathrm{mga}}{(a+b)}---\left(4\right)$

式(4)中，g表示重力加速度，m表示车辆质量，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离；

整车纵向动力学模型表述为：

F_x=F_xf+F_xr=ma_x+D_av²+C_rollmg  (5)

式(5)中，a_x表示车辆实时纵向加速度，由加速度计测量获得，D_a表示空气阻力常数，C_roll表示滚动阻力系数，v表示车辆绝对速度，F_x表示整车纵向力；

本发明采用简化魔术公式轮胎模型，下角标j取f或r，f或r分别表示前或后轮轴：

式(6)中，μ表示实时道路附着系数，B、C表示模型系数，假定各车轮的路面附着系数状况相同，前、后轮轴的归一化牵引力分别为：

由式(5)，(7)，(8)式可得

F_x=F_xf+F_xr=μ{F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}  (9)

3)基于最小二乘法的纵向道路附着系数初步估计

将式(9)表示为参数识别标准形式：

式(10)中，k表示离散时刻，y(k)=F_x表示***输出量，可由(5)式计算得知；θ(k)=μ_RLS表示待估参数矢量，其中μ_RLS表示由RLS方法所初步估计的道路附着系数；表示输入回归矢量，本发明中上角标′表示对矩阵转置；e(k)表示识别误差；则利用带遗忘因子的RLS算法实时确定唯一的未知量——道路附着系数的估计步骤如下：

(1)计算***输出变量y(k)，并计算迭代矢量

(2)计算识别误差e(k)；

(3)计算增益矢量N(k)；

其中，

其中，参数λ为遗忘因子；

(4)计算待估参数矢量θ(k)；

θ(k)=θ(k-1)+N(k)e(k)

至此，可实时初步估计出道路附着系数μ_RLS；

4)基于扩展卡尔曼滤波的纵向道路附着系数估计

将RLS方法所估计出的μ_RLS作为观测量，B、C和μ作为扩充状态，利用纵向动力学方程(5)建立扩展卡尔曼滤波模型，离散化后的状态方程与观测方程如式(11)：

X(k)=f(X(k-1))+W(k-1)

                  (11)

Z(k)=h(X(k))+V(k)

式(11)中，***状态向量X=[x₁，x₂，x₃，x₄]′，且x₁=v，x₂=μ，x₃=B，x₄=C；W表示零均值的***高斯白噪声向量，W=[w₁ w₂ w₃ w₄]′，其中w₁、w₂、w₃及w₄分别表示四个***高斯白噪声分量，W对应的***噪声协方差阵Q(k-1)为： $Q(k-1)=\left[\begin{array}{cccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2& 0\\ 0& 0& 0& {{\mathrm{\σ}}_{w_4}}^2\end{array}\right],$ 其中及分别表示***高斯白噪声w₁、w₂、w₃及w₄对应的方差；观测向量Z=[a_x-m，v_m，μ_RLS]′，a_x-m表示加速度传感器测量获得的车辆纵向加速度；v_m表示通过GPS测量获得的车辆速度；μ_RLS表示由递归最小二乘估计所得的道路附着系数初步估计值；V表示与W互不相关的零均值观测白噪声向量， 表示通过加速度传感器测量获得的车辆纵向加速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过GPS测量获得的车辆速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过上述最小二乘法估计获得的道路附着系数的观测噪声且是均值为0、方差为的高斯白噪声；V对应的观测噪声方差阵R可表示为 $R=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{a_{x-m}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{v_m}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{{\mathrm{\μ}}_{\mathrm{RLS}}}}^2\end{array}\right];$ f(·)和h(·)表示非线性的***函数向量和观测函数向量：

$f\left(X(k-1)\right)=\left[\begin{array}{c}{f}_1\left(X(k-1)\right)\\ f_2\left(X(k-1)\right)\\ f_3\left(X(k-1)\right)\\ f_4\left(X(k-1)\right)\end{array}\right],$ $h\left(X\left(k\right)\right)=\left[\begin{array}{c}{h}_1\left(X\left(k\right)\right)\\ h_2\left(X\left(k\right)\right)\\ h_3\left(X\left(k\right)\right)\end{array}\right],$ 其中

$f_1\left(X(k-1)\right)=v(k-1)-\frac{T\{D_a{\left[v(k-1)\right]}^2+C_{\mathrm{roll}}\mathrm{mg}\}}{m}+$

$\frac{\mathrm{T\μ}(k-1)\left\{F_{\mathrm{zf}}\sin \right[C(k-1)\arctan \left(B(k-1)s_f\right)]+F_{\mathrm{zr}}\sin [C(k-1)\arctan \left(B(k-1)S_r\right)\left]\right\}}{m}$

$,f_2\left(X(k-1)\right)=\mathrm{\μ}(k-1),f_3\left(X(k-1)\right)=B(k-1),f_4\left(X(k-1)\right)=C(k-1)$

$h_1\left(X\left(k\right)\right)=\frac{\mathrm{\μ}\left\{F_{\mathrm{zf}}\sin \right[C\arctan \left({\mathrm{Bs}}_f\right)]+F_{\mathrm{zr}}\sin [C\arctan \left({\mathrm{Bs}}_r\right)\left]\right\}-D_a{\left[v(k-1)\right]}^2-C_{\mathrm{roll}}\mathrm{mg}}{m}$

$,h_3\left(X\left(k\right)\right)=v,h_3\left(X\left(k\right)\right)=\mathrm{\μ}$

T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒；

对于式(11)所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，建立标准滤波递推过程，该递推过程包括时间更新和测量更新：

时间更新：

状态一步预测方程： $\hat{X}(k,k-1)=f\left(X(k-1)\right)$

一步预测误差方差阵P(k，k-1)：

P(k，k-1)=A(k，k-1)P(k-1)A′(k，k-1)+Q(k-1)

测量更新：

滤波增益矩阵K(k)：K(k)=P(k，k-1)·H′(k)·[H(k)P(k，k-1)H′(k)+R(k)]^-1

状态估计： $\hat{X}\left(k\right)=\hat{X}(k,k-1)+K\left(k\right)[Z\left(k\right)-h[\hat{X}(k,k-1)]]$

估计误差方差阵P(k)：P(k)=[I-K(k)·H(k)]·P(k，k-1)且I为4×4单位阵；

其中，A是***状态函数向量f对状态向量X求偏导数的雅可比矩阵，H是观测函数向量h对状态向量X求偏导数的雅可比矩阵，即矩阵A和H的第i行第j列元素A_[i，j](i＝1，2，3，4 j＝1，2，3，4)和H_[i，j](i＝1，2，3 j＝1，2，3，4)可分别通过下面的式子求得：

$A_{[i,j]}=\frac{{\∂f}_i}{{\∂x}_j}\left(\hat{X}(k,k-1)\right),(i=\mathrm{1,2,3,4},j=\mathrm{1,2,3,4})$

$H_{[i,j]}=\frac{{\∂h}_i}{{\∂x}_j}\left(\hat{X}(k,k-1)\right),(i=1,\mathrm{2,3},j=\mathrm{1,2,3,4})$

具体而言，各矩阵元素的取值如下：

A_[1，1]=1+T[(-2D_av)/m]

A_[1，2]=T{F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}/m

$A_{\left[\mathrm{1,3}\right]}=\mathrm{T\μ}[\frac{F_{\mathrm{zf}}{\mathrm{Cs}}_f\cos \left(C\arctan \left({\mathrm{Bs}}_f\right)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{\mathrm{zr}}{\mathrm{Cs}}_r\cos \left(C\arctan \left({\mathrm{Bs}}_r\right)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}]/m$

A_[1，4]=Tμ[F_zf(arctan(Bs_f)cos(Carctan(Bs_f)))+F_zr(arctan(Bs_r)cos(Carctan(Bs_r)))]/m

A_[2，1]=A_[2，3]=A_[2，4]=A_[3，1]=A_[3，2]=A_[3，4]=A_[4，1]=A_[4，2]=A_[4，3]=0

A_[2，2]=A_[3，3]A_[4，4]=1

H_[1，1]=(-2D_av)/m H_[1，2]={F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}/m

$H_{\left[\mathrm{1,3}\right]}=\mathrm{\μ}[\frac{F_{\mathrm{zf}}{\mathrm{Cs}}_f\cos \left(C\arctan \left({\mathrm{Bs}}_f\right)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{\mathrm{zr}}{\mathrm{Cs}}_r\cos \left(C\arctan \left({\mathrm{Bs}}_r\right)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}]/m$

H[_1，4]=μ[F_zf(arctan(Bs_f)cos(Carctan(Bs_f)))+F_zr(arctan(Bs_r)cos(Carctan(Bs_r)))]/m

H_[2，2]＝H_[2，3]＝H_[2，4]=H_[3，1]=H_[3，3]=H_[3，4]=0

H_[2，1]=H_[3，2]＝1

以上述EKF滤波递推所输出的滤波后的μ值作为最终道路附着系数估计结果。

本发明的优点及显著效果：

1.本发明提出了一种道路纵向附着系数估计方法，可用于平坦高速公路路面上汽车纵向主动安全控制对道路纵向附着系数的测量与估计需要，具有精度高、成本低、实时性和适应性好等优点。

2.本发明的方法是针对传统方法无法同时适应高滑移率和低滑移率工况，在整车纵向动力学模型和非线性简化魔术公式轮胎模型的基础上提出的，无论在高滑移率还是低滑移率工况下都可以获得较为准确的道路附着系数信息。

3.本发明提出的道路纵向附着系数估计方法对于道路附着系数的突变具有良好的适应性，响应时间短，能够满足平坦高速公路路面上汽车纵向主动安全控制的要求。

附图说明

图1是本发明所提出方法流程框图；

图2是简化魔术公式轮胎模型归一化牵引力与滑移率的关系；

图3是低滑移率单一附着系数工况下滑移率曲线；

图4是与图3对应的工况下本发明方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

图5是低滑移率附着系数突变工况下滑移率曲线；

图6是与图5对应的工况下本发明方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

图7是高滑移率单一附着系数工况下滑移率曲线；

图8是与图7对应的工况下本发明方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

图9是与图7对应的工况下传统滑移斜率方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

图10是高滑移率附着系数突变工况下滑移率曲线；

图11是与图10对应的工况下本发明方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

图12是与图10对应的工况下传统滑移斜率方法输出的附着系数估计结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；

具体实施方式

实施实例1

随着社会经济的发展，道路交通安全问题日益突出，并已成为全球性难题。全世界每年因交通事故都会造成大量的人员伤亡和财产损失，世界各国都在努力降低交通事故的发生。近年来，汽车主动安全技术得到了迅速的发展。汽车主动安全技术能够防患于未然，主动避免事故的发生，已成为现代汽车最主要的发展方向之一。目前常见的主动安全技术主要包括汽车防抱死制动***(ABS)，车辆电子稳定程序(ESP)，牵引力控制***(TCS)，电控驱动防滑***(ASR)，四轮转向稳定控制***(4WS)等。这些***效果的优劣很大程度上取决于能否“道路自适应”，即如果能够实时估计出道路附着系数，***就可以根据当前路况调节控制策略，提高车辆安全。纵向道路附着系数是影响汽车制动性能的一个关键因素，是碰撞预警/主动避撞(CW/CA)、自适应巡航(ACC)、防抱死制动***(ABS)等汽车纵向安全辅助***中必不可少的一个重要参数，实时估计纵向道路附着系数是设计符合实际要求的纵向安全算法的重要技术前提，其精度直接关系汽车的行驶安全性与稳定性，即纵向主动安全控制***能否有效工作在很大程度上依赖于道路附着系数能否被实时、准确的测量或估计。

目前，在汽车主动安全领域，道路附着系数主要分为直接测量和间接估计两类，直接测量方法是利用光、声、图像、雷达等传感器直接检测路面，测量一些对路面附着系数影响较大的因素，并根据以往经验预测当前道路附着系数的大小，但这些方法都需要额外加装传感器，且传感器成本都较高，难以实现大规模的商业应用；其次需要进行大量的测试训练，识别精度很大程度上依赖于经验，难以准确估算没有测试和训练过的路况的附着系数。间接估计方法是通过对汽车的运行过程进行运动学或动力学建模，结合轮胎模型，将有关低成本的车载传感器(如轮速传感器、陀螺仪、加速度计以及GPS等)信息作为观测信息，进而利用适当的滤波估计算法实现对道路附着系数的估计。已有的间接方法包括基于车辆侧向动力学和基于纵向动力学的研究两种，其中，前者主要用于汽车转向运动，而对于汽车纵向安全辅助***，较为关注的是基于车辆纵向动力学的附着系数估算。不同附着系数的路面上具有不同的滑移率和归一化牵引力关系，这是利用车辆纵向动力学估算附着系数的基础。基于这种关系，就可以通过估计滑移斜率，利用分类算法来估计附着系数，这是当前国内外学者所使用最多的方法。该方法认为小滑移率范围内的归一化牵引力与滑移率成正比，其比例关系称为滑移斜率，而道路附着系数与滑移斜率成正比，因此，通过估计滑移斜率，就能估计出轮胎路面的附着系数。但这种方法只能应用于低滑移率工况，而在高滑移率工况下，归一化牵引力与滑移斜率成正比的假设不再成立，因此无法进行道路附着系数估计。同时在同一路面上滑移斜率随其他因素的影响也存在一定差异，给附着系数的识别带来了一定的困难。另一种常用的方法是采用幂函数、指数函数、对数函数等不同形式的非线性公式及其组合来实时根据测量数据拟合归一化牵引力与滑移率曲线，从而根据计算所拟合曲线的峰值估计道路附着系数，但这种方法需要在较大轮胎滑移率条件下才能使估算结果更加接近于真实值，同时，该方法对归一化牵引力与滑移率曲线的拟合是基于所假设的非线性轮胎模型的，若选用轮胎模型过于复杂，则实时性难以保证，若选取的轮胎模型过于简单，则又无法保证准确性。可见，当前在附着系数估计上所存在的主要问题在于难以同时在高滑移率和低滑移率工况下都取得较好的估计效果，限制了其在汽车纵向主动安全***上的应用。

为实现在平坦高速公路路面的高滑移率和低滑移率工况下都能对道路附着系数进行准确、可靠估计，本发明提出了一种道路纵向附着系数估计方法。这种方法基于基于递归最小二乘(Recursive least squares，RLS)和扩展卡尔曼滤波(Extended Kalman Filter，EKF)。本发明针对现有方法无法同时适用于高滑移率和低滑移率工况，采用更符合实际车辆运行过程的简化魔术公式轮胎模型，结合纵向整车动力学模型利用递归最小二乘法进行道路附着系数估计，同时充分利用低成本车载传感器信息作为***的观测量，通过扩展卡尔曼滤波方法实现轮胎模型系数的自适应调整，并去除纵向道路附着系数估计结果所含噪声，输出最终道路纵向附着系数估计结果，具有精度高、成本低、实时性和适应性好等显著优点，所估计出的道路附着系数可用于平坦高速公路路面上行驶的汽车纵向主动安全控制。本发明的具体思路如下：

递归最小二乘是对未知矢量的迭代算法，以模型误差的最小方差为目标，对于每个采样周期，使用已有采样数据通过反复迭代计算未知矢量。卡尔曼滤波器是以最小均方差为准则的最优状态估计滤波器，它不需要储存过去的测量值，只根据当前的观测值和前一时刻的估计值，利用计算机进行递推计算，便可实现对实时信号的估计。递归最小二乘和卡尔曼滤波都具有数据存储量小、算法简便的特点。

为适应在平坦高速公路路面的高滑移率和低滑移率工况下汽车纵向主动安全控制对纵向道路附着系数的测量与估计要求，首先对汽车和轮胎进行适当的动力学建模。针对本发明的应用领域，本发明对于行驶在通常道路交通环境上的前轮转向的四轮车辆(目前应有最广的情况，典型例子如前轮转向的轿车)，忽略侧向、横摆、俯仰、侧倾和上下弹跳运动，忽略左右轮差异，由于本发明的范围为平坦路面，故可忽略道路坡度，假定各车轮的路面附着系数状况相同，整车纵向动力学模型表述为：

ma_x=F_x-D_av²-C_rollmg     (1)

式(1)中，m表示车辆总质量，a_x表示车辆实时纵向加速度，由加速度传感器获取，F_x表示整车纵向力，D_a表示空气阻力常数，C_roll表示滚动阻力系数，v表示车辆绝对速度，由GPS获取，g表示重力加速度。

式(1)可变形为：

F_x=F_xf+F_xr=ma_x+D_av²+C_rollmg     (2)

式(2)中，F_xf和F_xr分别表示前、后轮轴上纵向力，F_xf和F_xr可统一记为F_xj(j＝f，r)。

魔术公式轮胎模型是公认的拟合精度最高的经验轮胎模型，无论是高滑移率还是低滑移率工况都同样适用，但它是由三角函数组合而成的复杂的非线性函数，且模型中需要确定的未知因子也较多，计算量较大，不适宜于实时使用。因此，本文采用其简化模型[可参考文献：边明远，用于纵向道路附着系数评估的简化轮胎模型[J].重庆理工大学学报(自然科学)，2012,26(1)：1-5.]。用表示归一化牵引力，即又可分为前轮轴归一化牵引力和后轮轴归一化牵引力下角标j取f或r，f或r分别表示前或后轮轴：

(3)

式(3)中，F_xf和F_xr分别表示前、后轮轴上纵向力，F_xf和F_xr可统一记为F_xj(j＝f，r)，F_zf和F_zr分别表示分配到前或后轮轴上的垂向载荷，可统一记为F_zj(j＝f，r)且可按下式计算：

$\begin{array}{cc}{F}_{\mathrm{zf}}=\frac{\mathrm{mgb}}{(a+b)},& F_{\mathrm{zr}}=\frac{\mathrm{mga}}{(a+b)}\end{array}$

(4)

式(4)中，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离。

用s_j(j＝f，r)表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_f和后轮轴纵向滑移率s_r，下角标j取f或r，f或r分别表示前或后轮轴，s_j计算方法为：

(5)

式(5)中，R_e表示车轮轮胎半径；ν_tf和ν_tr分别表示前、后轮轴上沿轮胎方向的速度，v_tf和v_tr可统一记为v_tj(j＝f，r)，由于本发明所适用的场景为高速公路，高速公路的道路曲率较小，车辆行驶于高速公路上时，其横摆角速度也较小，且本发明所用的方法是基于纵向动力学以进行纵向道路附着系数估计，因此本发明中v_tj可近似等于车辆绝对速度v，通过GPS测量获取；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_j(j＝f，r)且

${\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{fR}}+{\mathrm{\ω}}_{\mathrm{fL}})$

                     (6)

${\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{rR}}+{\mathrm{\ω}}_{\mathrm{rL}})$

式(6)中，ω_fL、ω_fR、ω_rL和ω_rR分别表示左前轮、右前轮、左后轮和右后轮的旋转角速度，通过利用四个轮速传感器测量获得；

简化魔术公式轮胎模型如图2所示，图2中为μ＝0.2，0.4，0.6，0.8，1.0时归一化牵引力与滑移率的关系，可见在低滑移率时成近似线性关系，这也是传统滑移斜率方法利用滑移斜率估计道路附着系数的基础。根据简化魔术公式轮胎模型，有

式(7)中，μ表示实时道路附着系数，B、C表示模型系数，前、后轮轴的归一化牵引力分别为：

由式(2)，(8)，(9)式可得

F_x=F_xf+F_xr=μ{F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}   (10)

基于最小二乘法的纵向道路附着系数初步估计步骤如下：

将式(10)表示为参数识别标准形式：

式(11)中，k表示离散时刻，y(k)=F_x表示***输出量，可由(2)式计算得知，θ(k)=μ_RLS表示待估参数矢量，其中μ_RLS表示由RLS方法所初步估计的道路附着系数；表示输入回归矢量，本发明中上角标′表示对矩阵转置；e(k)表示识别误差，则利用带遗忘因子的RLS算法实时确定唯一的未知量——道路附着系数的估计步骤如下：

(1)计算***输出变量y(k)，并计算迭代矢量

(2)计算识别误差e(k)：

(3)计算增益矢量N(k)：

其中，

其中，参数λ为遗忘因子，能有效减少不再与模型相关的旧数据的影响，并防止协方差发散，通常取值范围在[0.9，1]，本发明取0.995。

(4)计算待估参数矢量θ(k)

θ(k)＝θ(k-1)+N(k)e(k)

至此，可实时初步估计出道路附着系数μ_RLS。

在车辆运行过程中，轮胎模型参数B和C值并非一直为常量，受车速等因素影响会有一定改变，同时，所测得的整车纵向车速也含有较大的噪声，都会对附着系数的估计产生一定的影响。基于这样的考虑，将RLS所估计出的μ_RLS作为观测量，B、C和μ作为扩充状态，利用纵向动力学方程(1)建立扩展卡尔曼滤波模型，实现B、C值的实时更新和v、μ的滤波，离散化后的状态方程与观测方程如式(12)：

X(k)＝f(X(k-1))+W(k-1)

                    (12)

Z(k)＝h(X(k))+V(k)

式(12)中，***状态向量X＝[x₁，x₂，x₃，x₄]′，且x₁＝v，x₂＝μ，x₃＝B，x₄＝C；本发明中上角标′表示对矩阵转置；W表示零均值的***高斯白噪声向量，W＝[w₁ w₂ w₃ w₄]′，其中w₁、w₂、w₃及w₄分别表示四个***高斯白噪声分量，W对应的***噪声协方差阵Q(k-1)为： $Q(k-1)=\left[\begin{array}{cccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2& 0\\ 0& 0& 0& {{\mathrm{\σ}}_{w_4}}^2\end{array}\right],$ 其中及分别表示***高斯白噪声w₁、w₂、w₃及w₄对应的方差；观测向量Z＝[a_x-m，v_m，μ_RLS]′，a_x-m表示加速度传感器测量获得的车辆纵向加速度；v_m表示通过GPS测量获得的车辆速度；μ_RLS表示由递归最小二乘估计所得的道路附着系数初步估计值；V表示与W互不相关的零均值观测白噪声向里， 表示通过加速度传感器测量获得的车辆纵向加速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过GPS测量获得的车辆速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过上述最小二乘法估计获得的道路附着系数的观测噪声且是均值为0、方差为的高斯白噪声；V对应的观测噪声方差阵R可表示为 $R=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{a_{x-m}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{v_m}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{{\mathrm{\μ}}_{\mathrm{RLS}}}}^2\end{array}\right];$ f(·)和h(·)表示非线性的***函数向量和观测函数向量：

$f\left(X(k-1)\right)=\left[\begin{array}{c}{f}_1\left(x(k-1)\right)\\ f_2\left(x(k-1)\right)\\ f_3\left(x(k-1)\right)\\ f_4\left(x(k-1)\right)\end{array}\right],$ $h\left(X\left(k\right)\right)=\left[\begin{array}{c}{h}_1\left(X\left(k\right)\right)\\ h_2\left(X\left(k\right)\right)\\ h_3\left(X\left(k\right)\right)\end{array}\right],$ 其中

$f_1\left(X(k-1)\right)=v(k-1)-\frac{T\{D_a{\left[v(k-1)\right]}^2+C_{\mathrm{roll}}\mathrm{mg}\}}{m}+$

$\frac{\mathrm{T\μ}(k-1)\left\{F_{\mathrm{zf}}\sin \right[C(k-1)\arctan \left(B(k-1)s_f\right)]+F_{\mathrm{zr}}\sin [C(k-1)\arctan \left(B(k-1)s_r\right)\left]\right\}}{m}$

$,f_2\left(X(k-1)\right)=\mathrm{\μ}(k-1),f_3\left(X(k-1)\right)=B(k-1),f_4\left(X(k-1)\right)=C(k-1)$

$h_1\left(X\left(k\right)\right)=\frac{\mathrm{\μ}\left\{F_{\mathrm{zf}}\sin \right[C\arctan \left({\mathrm{Bs}}_f\right)]+F_{\mathrm{zr}}\sin [C\arctan \left({\mathrm{Bs}}_r\right)\left]\right\}-D_a{\left[v(k-1)\right]}^2-C_{\mathrm{roll}}\mathrm{mg}}{m}$

$,h_2\left(X\left(k\right)\right)=v,h_3\left(X\left(k\right)\right)=\mathrm{\μ}$

T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒。

对于式(12)所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，建立标准滤波递推过程，该递推过程包括时间更新和测量更新：

时间更新：

状态一步预测方程： $\hat{X}(k,k-1)=f\left(X(k-1)\right)$

一步预测误差方差阵P(k，k-1)：

P(k，k-1)＝A(k，k-1)P(k-1)A′(k，k-1)+Q(k-1)

测量更新：

滤波增益矩阵K(k)：K(k)＝P(k，k-1)·H′(k)·[H(k)P(k，k-1)H′(k)+R(k)]^-1

状态估计： $\hat{X}\left(k\right)=\hat{X}(k,k-1)+K\left(k\right)[Z\left(k\right)-h[\hat{X}(k,k-1)\left]\right]$

估计误差方差阵P(k)：P(k)＝[I-K(k)·H(k)]·P(k，k-1)且I为4×4单位阵；

其中，A是***状态函数向量f对状态向量X求偏导数的雅可比矩阵，H是观测函数向量h对状态向量X求偏导数的雅可比矩阵，即矩阵A和H的第i行第j列元素A_[i，j](i＝1，2，3，4 j＝1，2，3，4)和H_[i，j](i＝1，2，3 j＝1，2，3，4)可分别通过下面的式子求得：

$A_{[i,j]}=\frac{{\∂f}_i}{{\∂x}_j}\left(\hat{X}(k,k-1)\right),(i=\mathrm{1,2,3,4},j=\mathrm{1,2,3,4})$

$H_{[i,j]}=\frac{{\∂h}_i}{{\∂x}_j}\left(\hat{X}(k,k-1)\right),(i=\mathrm{1,2,3},j=\mathrm{1,2,3,4})$

具体而言，各矩阵元素的取值如下：

A_[1，1]＝1+T[(-2D_av)/m]

A_[1，2]＝T{F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}/m

$A_{\left[\mathrm{1,3}\right]}=\mathrm{T\μ}[\frac{F_{\mathrm{zf}}{\mathrm{Cs}}_f\cos \left(C\arctan \left({\mathrm{Bs}}_f\right)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{\mathrm{zr}}{\mathrm{Cs}}_r\cos \left(C\arctan \left({\mathrm{Bs}}_r\right)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}]/m$

A_[1，4]＝Tμ[F_zf(arctan(Bs_f)cos(Carctan(Bs_f)))+F_zr(arctan(Bs_r)cos(Carctan(Bs_r)))]/m

A_[2，1]＝A_[2，3]＝A_[2，4]＝A_[3，1]＝A_[3，2]＝A_[3，4]＝A_[4，1]＝A_[4，2]＝A_[4，3]＝0

A_[2，2]＝A_[3，3]＝A_[4，4]＝1

H_[1，1]＝(-2D_av)/m H_[1，2]＝{F_zfsin[Carctan(Bs_f)]+F_zrsin[Carctan(Bs_r)]}/m

$H_{\left[\mathrm{1,3}\right]}=\mathrm{\μ}[\frac{F_{\mathrm{zf}}{\mathrm{Cs}}_f\cos \left(C\arctan \left({\mathrm{Bs}}_f\right)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{\mathrm{zr}}{\mathrm{Cs}}_r\cos \left(C\arctan \left({\mathrm{Bs}}_r\right)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}]/m$

H_[1，4]＝μ[F_zf(arctan(Bs_f)cos(Carctan(Bs_f)))+F_zr(arctan(Bs_r)cos(Carctan(Bs_r)))]/m

H_[2，2]＝H_[2，3]＝H_[2，4]＝H_[3，1]＝H_[3，3]＝H_[3，4]＝0

H_[2，1]＝H_[3，2]＝1

在上述滤波递推计算过程中，可确定汽车在每个时刻的轮胎模型参数B、C、速度v和道路附着系数μ，将EKF所输出的状态量B、C反馈回轮胎模型，形成闭环反馈，实时更新递归最小二乘估计，可进一步提高附着系数估计精度，及时响应行驶状况的改变，以EKF所输出的道路附着系数μ值即作为最终估计结果。

实施实例2

为检验本发明提出的基于改进的扩展卡尔曼滤波的车辆运行状态估计方法的实际效果，在专业的汽车动力学仿真软件CarSim上进行了仿真验证实验。

CarSim是由美国MSC(Mechanical Simulation Corporation)公司开发的专门针对车辆动力学的仿真软件，目前已被国际上众多的汽车制造商、零部件供应商所采用，被广泛地应用于现代汽车控制***的商业开发，已成为汽车行业的标准软件，享有很高的声誉。Carsim内的车辆动力学模型是通过分别对汽车的车体、悬架、转向、制动等各子***以及各个轮胎的高逼真建模来实现的，具有很高的自由度，能够提供非常接近实际的准确的车辆运行状态信息，因此，Carsim输出的车辆运行状态信息可作为车辆的参考输出。

为检验本发明提出的算法在低滑移率和高滑移率工况下的估计效果，以及本发明提出的算法对道路附着系数突变的响应，仿真实验中分别设置了高滑移率和低滑移率工况，以及附着系数的突变，并在高滑移率工况下与传统的滑移斜率方法进行了对比。基于滑移斜率的方法认为小滑移率范围内的归一化牵引力与滑移率成正比，其比例关系称为滑移斜率，而道路附着系数与滑移斜率成正比，因此，通过估计滑移斜率，就能估计出轮胎路面的附着系数。在低滑移率工况下，仿真时间设置为45s，在高滑移率工况下，仿真时间设置为50s。采样频率为10Hz，即周期T为0.1s。所用车辆是一个前轮转向的四轮车，主要参数如下：m＝1220(千克)、R_e＝310.8(毫米)、a＝1.040(米)、b＝1.560(米)。设定四个车轮的线速度(通过轮速传感器测得的角速度乘以轮胎半径得到)的测量噪声均为均值是0、标准差是0.04(米/秒)的高斯白噪声。轮胎模型参数B、C分别取初值为14和1.3[可参考文献：Fredrik Gustafsson，Automotive Safety Systems.Replacing costly sensorswith software algorithms，2009，IEEE Signal Processing Magazine，2009，26(4)：32-47.]，重力加速度g取9.78，递归最小二乘的遗忘因子设置为0.995，卡尔曼滤波的***零均值高斯白噪声的标准差分别为(米/秒)、及卡尔曼滤波的三个观测量的零均值高斯白噪声的标准差分别为(米/秒²)、(米/秒)及有关结果如图3～图12所示。

图3-图6为低滑移率工况下仿真结果，其中图3为附着系数设置为0.2时的滑移率，图4为图3对应工况下利用本发明算法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)，图5为附着系数设置为从0.8跃变至0.2的滑移率，图6为图5对应工况下利用本发明算法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)，由图3、图5可见其滑移率都较小，可见本文所提出的方法在低滑移率工况下，能够迅速识别道路附着系数，准确性较高，误差小于0.1，由图6可见，在道路附着系数跃变时，算法能在2s内迅速收敛到真值附近，快速响应道路附着系数的突变。

通过设置车辆采取较为强烈的制动操作以产生高滑移率工况。并以传统的基于滑移斜率的方法估计道路附着系数以进行比较，基于滑移斜率的方法认为小滑移率范围内的归一化牵引力与滑移率成正比，其比例关系称为滑移斜率，而道路附着系数与滑移斜率成正比，因此，通过估计滑移斜率，就能估计出轮胎路面的附着系数。图7-图12为在高滑移率工况下的仿真结果。其中图7为附着系数设置为0.2时的滑移率，图8为图7对应工况下利用本发明算法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)，图9为图7对应工况下利用传统滑移斜率方法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)；图10为附着系数设置为从0.8跃变至0.2的滑移率，图11为图10对应工况下利用本发明算法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)，图12为图10对应工况下利用传统滑移斜率方法进行附着系数估计的结果(图中附着系数估计结果用黑实线代表，Carsim输出真值用黑虚线代表)，由图7、图10可见纵向滑移率都较大，由于纵向滑移率不再近似与轮胎纵向力成线性关系，传统滑移斜率方法产生了极大的误差，而本发明所采用方法在在估计的准确性、对路面突变响应的及时性方面具有很好的效果。

Claims (1)

1.一种道路纵向附着系数估计方法，其特征在于：在平坦的高速公路路面上，针对前轮转向四轮汽车，基于整车纵向动力学模型和简化魔术公式轮胎模型，利用带遗忘因子的递归最小二乘(Recursive least squares,RLS)方法实时初步估计出道路纵向附着系数，进一步将初步估计出的道路纵向附着系数与轮胎模型参数作为扩充状态，利用扩展卡尔曼滤波(Extended Kalman Filter,EKF)算法，滤除信号噪声，并实现轮胎模型参数的自适应调整，最终实时获取准确、鲁棒的道路纵向附着系数估计，能够同时适应平坦高速公路上的高滑移率和低滑移率工况；

具体步骤包括：

1)计算纵向滑移率

用s_j(j＝f,r)表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_f和后轮轴纵向滑移率s_r，下角标j取f或r，f、r分别表示前、后轮轴，s_j计算方法为：

式(1)中，R_e表示轮胎半径；v_tf和v_tr分别表示前、后轮轴上沿轮胎方向的速度，以车辆绝对速度v代替，通过GPS测量获取；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_j(j＝f,r)且

${\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{fR}+{\mathrm{\ω}}_{fL})$                    (2)

${\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{rR}+{\mathrm{\ω}}_{rL})$

式(2)中，ω_fL、ω_fR、ω_rL和ω_rR分别表示左前轮、右前轮、左后轮和右后轮的旋转角速度，通过四个轮速传感器测量获得；

2)计算整车纵向力

用表示归一化牵引力，即又可分为前轮轴归一化牵引力和后轮轴归一化牵引力下角标j取f或r，f、r分别表示前、后轮轴：

式(3)中，F_xf、F_xr分别表示前、后轮轴上纵向力，F_xf和F_xr可统一记为F_xj(j＝f,r)，F_zf、F_zr分别表示分配到前、后轮轴上的垂向载荷，可统一记为F_zj(j＝f,r)且可按下式计算：

$F_{zf}=\frac{mgb}{(a+b)},F_{zr}=\frac{mga}{(a+b)}---\left(4\right)$

式(4)中，g表示重力加速度，m表示车辆质量，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离；

整车纵向动力学模型表述为：

F_x＝F_xf+F_xr＝ma_x+D_av²+C_rollmg             (5)

式(5)中，a_x表示车辆实时纵向加速度，由加速度计测量获得，D_a表示空气阻力常数，C_roll表示滚动阻力系数，v表示车辆绝对速度，F_x表示整车纵向力；

本发明采用简化魔术公式轮胎模型，下角标j取f或r，f、r分别表示前、后轮轴：

式(6)中，μ表示实时道路纵向附着系数，B、C表示轮胎模型参数轮胎模型参数B、C分别取初值为14和1.3，假定各车轮的路面附着系数状况相同，前、后轮轴的归一化牵引力分别为：

由式(5)，(7)，(8)式可得

F_x＝F_xf+F_xr＝μ{F_zf sin[C arctan(Bs_f)]+F_zr sin[C arctan(Bs_r)]}       (9)

3)基于最小二乘法的道路纵向附着系数初步估计

将式(9)表示为参数识别标准形式：

式(10)中，k表示离散时刻，y(k)＝F_x表示***输出量，可由(5)式计算得知；θ(k)＝μ_RLS表示待估参数矢量，其中μ_RLS表示由递归最小二乘方法所初步估计的道路纵向附着系数；表示输入回归矢量，上角标'表示对矩阵转置；e(k)表示识别误差；则利用带遗忘因子的递归最小二乘算法实时确定唯一的未知量——道路纵向附着系数的估计步骤如下：

(1)计算***输出变量y(k)，并计算迭代矢量

(2)计算识别误差e(k)；

(3)计算增益矢量N(k)；

其中，

其中，参数λ为遗忘因子；

(4)计算待估参数矢量θ(k)；

θ(k)＝θ(k-1)+N(k)e(k)

至此，可实时初步估计出道路纵向附着系数μ_RLS；

4)基于扩展卡尔曼滤波的道路纵向附着系数估计

将递归最小二乘方法所初步估计出的μ_RLS作为观测量，B、C和μ作为扩充状态，利用纵向动力学方程(5)建立扩展卡尔曼滤波模型，离散化后的状态方程与观测方程如式(11)：

X(k)＝f(X(k-1))+W(k-1)                 (11)

Z(k)＝h(X(k))+V(k)

式(11)中，***状态向量X＝[x₁,x₂,x₃,x₄]′，且x₁＝v，x₂＝μ，x₃＝B,x₄＝C；W表示零均值的***高斯白噪声向量，W＝[w₁,w₂,w₃,w₄]′，其中w₁、w₂、w₃及w₄分别表示四个***高斯白噪声分量，W对应的***噪声协方差阵Q(k-1)为： $Q(k-1)=\left[\begin{array}{cccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2& 0\\ 0& 0& 0& {{\mathrm{\σ}}_{w_4}}^2\end{array}\right],$ 其中及分别表示***高斯白噪声w₁、w₂、w₃及w₄对应的方差；观测向量Z＝[a_x,v,μ_RLS]′，μ_RLS表示由递归最小二乘估计所得的道路纵向附着系数初步估计值；V表示与W互不相关的零均值观测白噪声向量， 表示通过加速度传感器测量获得的车辆纵向加速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过GPS测量获得的车辆速度的观测噪声且是均值为0、方差为的高斯白噪声；表示通过上述最小二乘法估计获得的道路纵向附着系数的观测噪声且是均值为0、方差为的高斯白噪声；V对应的观测噪声方差阵R可表示为 $R=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{a_{x-m}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{v_m}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{{\mathrm{\μ}}_{RLS}}}^2\end{array}\right];$ f(·)和h(·)表示非线性的***函数向量和观测函数向量：

$f\left(X\right(k-1\left)\right)=\left[\begin{array}{c}f\left(X\left(k-1\right)\right)\\ f_2\left(X\left(k-1\right)\right)\\ f_3\left(X\left(k-1\right)\right)\\ f_4\left(X\left(k-1\right)\right)\end{array}\right],h\left(X\right(k\left)\right)=\left[\begin{array}{c}h\left(X\left(k\right)\right)\\ h_2\left(X\left(k\right)\right)\\ h_3\left(X\left(k\right)\right)\end{array}\right],$ 其中

$\begin{array}{c}{f}_1\left(X\right(k-1\left)\right)=v\left(k-1\right)-\frac{T\left\{D_a{\[v\left(k-1\right)\]}^2+C_{roll}mg\right\}}{m}+\\ \frac{T\mathrm{\μ}(k-1)\{F_{zf}\sin \[C(k-1)arctan\left(B\right(k-1\left)s_f\right)\]+F_{zr}sin\[C(k-1)arctan\left(B\right(k-1\left)s_r\right)\]\}}{m}\\ ,f_2\left(X\right(k-1\left)\right)=\mathrm{\μ}(k-1),f_3\left(X\right(k-1\left)\right)=B(k-1),f_4\left(X\right(k-1\left)\right)=C(k-1)\end{array}$

$\begin{array}{c}{h}_1\left(X\right(k\left)\right)=\frac{\mathrm{\μ}\{F_{zf}sin\[Carctan\left({\mathrm{Bs}}_f\right)\]+F_{zr}sin\[Carctan\left({\mathrm{Bs}}_r\right)\]\}-D_a{v}^2-C_{roll}mg}{m},h_2\left(X\right(k\left)\right)=v\\ h_3\left(X\right(k\left)\right)=\mathrm{\μ}\end{array},$

T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒；

对于式(11)所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，建立标准滤波递推过程，该递推过程包括时间更新和测量更新：

时间更新：

状态一步预测方程： $\hat{X}(k,k-1)=f\left(X\right(k-1\left)\right)$

一步预测误差方差阵P(k,k-1)：

P(k,k-1)＝A(k,k-1)P(k-1)A′(k,k-1)+Q(k-1)

测量更新：

滤波增益矩阵K(k)：K(k)＝P(k,k-1)H′(k)[H(k)P(k,k-1)H′(k)+R(k)]^-1

状态估计： $\hat{X}\left(k\right)=\hat{X}(k,k-1)+K\left(k\right)\[Z\left(k\right)-h\[\hat{X}(k,k-1)\]\]$

估计误差方差阵P(k)：P(k)＝[I-K(k)H(k)]P(k,k-1)且I为4×4单位阵；

其中，A是***状态函数向量f对状态向量X求偏导数的雅可比矩阵，H是观测函数向量h对状态向量X求偏导数的雅可比矩阵，即矩阵A和H的第i行第j列元素A_[i,j](i＝1,2,3,4j＝1,2,3,4)和H_[i,j](i＝1,2,3  j＝1,2,3,4)可分别通过下面的式子求得：

$A_{\[i,j\]}=\frac{\∂f_i}{\∂x_j}\left(\hat{X}\right(k,k-1\left)\right),(i=1,2,3,4,j=1,2,3,4)$

$H_{\[i,j\]}=\frac{\∂h_i}{\∂x_j}\left(\hat{X}\right(k,k-1\left)\right),(i=1,2,3,j=1,2,3,4)$

具体而言，各矩阵元素的取值如下：

A_[1,1]＝1+T[(-2D_av)/m]

A_[1,2]＝T{F_zf sin[C arctan(Bs_f)]+F_zr sin[C arctan(Bs_r)]}/m

$A_{\[1,3\]}=T\mathrm{\μ}\[\frac{F_{zf}{\mathrm{Cs}}_f\cos \left(C\arctan \right({\mathrm{Bs}}_f\left)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{zr}{\mathrm{Cs}}_r\cos \left(C\arctan \right({\mathrm{Bs}}_r\left)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}\]/m$

A_[1,4]＝Tμ[F_zf(arctan(Bs_f)cos(C arctan(Bs_f)))+F_zr(arctan(Bs_r)cos(C arctan(Bs_r)))]/m

A_[2,1]＝A_[2,3]＝A_[2,4]＝A_[3,1]＝A_[3,2]＝A_[3,4]＝A_[4,1]＝A_[4,2]＝A_[4,3]＝0

A_[2,2]＝A_[3,3]＝A_[4,4]＝1

H_[1,1]＝(-2D_av)/m H_[1,2]＝{F_zf sin[C arctan(Bs_f)]+F_zr sin[C arctan(Bs_r)]}/m

$H_{\[1,3\]}=\mathrm{\μ}\[\frac{F_{zf}{\mathrm{Cs}}_f\cos \left(C\arctan \right({\mathrm{Bs}}_f\left)\right)}{1+{\left({\mathrm{Bs}}_f\right)}^2}+\frac{F_{zr}{\mathrm{Cs}}_r\cos \left(C\arctan \right({\mathrm{Bs}}_r\left)\right)}{1+{\left({\mathrm{Bs}}_r\right)}^2}\]/m$

H_[1,4]＝μ[F_zf(arctan(Bs_f)cos(C arctan(Bs_f)))+F_zr(arctan(Bs_r)cos(C arctan(Bs_r)))]/m

H_[2,2]＝H_[2,3]＝H_[2,4]＝H_[3,1]＝H_[3,3]＝H_[3,4]＝0

H_[2,1]＝H_[3,2]＝1

以上述扩展卡尔曼滤波递推所输出的滤波后的μ值作为最终道路纵向附着系数估计结果。