CN103434511A - 一种车速与道路附着系数的联合估计方法 - Google Patents

Abstract

本发明公开了一种车速与道路附着系数的联合估计方法。本方法基于非线性整车动力学模型和轮胎纵向力模型，在不同道路附着系数条件下，分别建立不同的多个卡尔曼滤波模型，同时利用车载轮速和方向盘转角传感器信息来确定各卡尔曼滤波系统的外部输入量和观测量，进一步通过交互多模型算法实现不同滤波系统的交互，从而实现不同道路附着系数条件下对车辆纵向、侧向车速的自适应估计，并根据交互多模型算法中计算出的各卡尔曼滤波模型的模型概率实现道路附着系数的实时估计，达到全面自适应的效果。

Description

一种车速与道路附着系数的联合估计方法

技术领域

本发明涉及一种车速与道路附着系数的联合估计方法，其目的在于实现汽车运行过程中对于车辆纵向速度、侧向速度以及道路附着系数的实时估计，这些估计值可用于汽车主动安全的相关控制，属于汽车主动安全测量及控制领域。

背景技术

随着社会经济的发展，道路交通安全问题日益突出，并已成为全球性难题。全世界每年因交通事故都会造成大量的人员伤亡和财产损失，世界各国都在努力降低交通事故的发生。近年来，汽车主动安全技术得到了迅速的发展。汽车主动安全技术能够防患于未然，主动避免事故的发生，已成为现代汽车最主要的发展方向之一。目前常见的主动安全技术主要包括防抱死制动系统（ABS），车辆电子稳定程序（ESP），牵引力控制系统（TCS），电控驱动防滑系统（ASR），四轮转向稳定控制系统（4WS）等。这些系统通常涉及汽车的纵向前进速度、侧向速度以及质心侧偏角等运行状态的测量或估计，这些运行状态的测量或估计可用于后续的汽车主动安全控制，因此其精度直接关系汽车的行驶安全性与稳定性。追尾碰撞预警/避免系统（CW/CA）、制动防抱死系统（ABS）等纵向主动安全系统依赖于纵向车速的准确估计，而以电子稳定程序（ESP）、四轮转向控制（4WS）为代表的侧向安全系统则依赖于对质心侧偏角的准确估计，而质心侧偏角可根据纵向车速和侧向车速计算得知。所以说，纵向和侧向车速信息获取的准确性和可靠性，直接决定着这些主动安全系统的有效性。同时，这些主动安全系统效果的优劣很大程度上取决于能否“道路自适应”，即如果能够实时估计出道路附着系数，系统就可以根据当前路况调节控制策略，提高车辆安全。因此，纵向、侧向车速与道路附着系数作为汽车主动安全的关键参数，对其进行实时、准确的测量或估计，是上述汽车主动安全系统发挥作用的重要前提与基础。

但是由于成本、技术等方面的原因，纵、侧向车速和道路附着系数信息通常难以直接测量，而是使用量产车上已装备的传感器，通过建立车辆运动或动力学模型，利用车辆状态估计技术来获取。

在车速估计方面，主要有基于运动学模型和基于动力学模型的估计方法，其中，基于运动学模型的估计方法虽然鲁棒性较好，不受模型参数影响，但对于传感器精度要求较高，车载低成本传感器的精度往往难以满足其需求，给其应用带来了很大的限制。而基于动力学模型的估计方法则对传感器精度要求相对较低，适宜于与车载传感器结合使用，因此成为了目前较为普遍的方法。在三自由度整车动力学模型或两自由度线性整车模型的基础上，结合线性或各种非线性轮胎模型，卡尔曼滤波算法、滑模观测器,模糊观测器,神经网络观测器及各种非线性观测器被广泛应用于车速估计中，取得了较好的效果。但基于动力学模型的方法对于模型本身的精度要求很高，要求模型尽可能准确的反应车辆的动力学特性，因而对模型参数准确性要求较高，同时对参数的变化较为敏感，道路附着系数作为重要的轮胎力学参数，在这些方法中往往被假设为已知且定常，从而根据事先设定的附着系数经验值，建立单一的轮胎模型计算轮胎力以估计车速，事实上，道路附着系数往往偏离事先预设的经验值，且在车辆运行过程中会变化而非定常，从而使模型不能很好的适应于当前道路情况，或导致车速估计结果无法及时响应道路条件的改变，严重影响其准确性。

道路附着系数一般无法直接测量，同样是通过车辆状态估计的方法来获取。现有的道路附着系数估计方法包括基于车辆侧向动力学和纵向动力学的方法两类，但这些估计方法或有其特定的适用范围（例如仅适宜于滑移率较低的线性区域），或需要一定的经验性，或存在难以测量的参数，或在实际使用中的效果有待进一步确认，同时，道路附着系数的准确估计也依赖于对于纵向、横向车速等车辆运行状态的准确估计，这些都限制了这些估计方法在汽车主动安全系统上的应用。

总而言之，由于道路附着系数和车辆纵向、横向车速在车辆运行过程中相互作用，相互影响，传统估计方法中对于道路附着系数定常的假设无法满足不同道路附着系数条件下对车速准确估计的要求，也无法实时估计出道路附着系数值。

发明内容

为实现在不同道路附着系数条件对车辆纵向速度和侧向速度的准确、可靠估计，并实时估计出相应的道路附着系数，本发明提出了一种基于交互多模型的车速与道路附着系数的联合估计方法。本发明提出的方法针对汽车在不同道路附着系数下的行驶工况，建立多个扩展卡尔曼滤波模型，同时充分利用低成本的车载轮速和方向盘转角传感器信息来建立滤波系统的外部输入量和观测量，进而通过交互多模型-扩展卡尔曼滤波算法，实现不同道路附着系数条件下对车辆纵向、侧向车速的自适应估计，并根据交互多模型算法中计算出的各卡尔曼滤波模型的模型概率实现道路附着系数的实时估计，具有精度高、成本低、实时性好、全面自适应等特点。

一种车速与道路附着系数的联合估计方法，其特征在于：本方法是针对前轮转向四轮汽车，基于非线性整车动力学模型和轮胎纵向力模型，在不同道路附着系数条件下，分别建立不同的多个卡尔曼滤波模型，同时利用车载轮速和方向盘转角传感器信息来确定建立各卡尔曼滤波系统的外部输入量和观测量。进一步通过交互多模型算法实现不同道路附着系数条件下对车辆纵向、侧向车速的自适应估计，并根据交互多模型算法中计算出的各卡尔曼滤波模型的模型概率实现道路附着系数的实时估计，达到全面自适应的效果；

具体步骤包括：

1）建立扩展卡尔曼滤波的状态方程和观测方程

针对道路附着系数分别为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0，即μ_j＝10×j（j=1,2,...10）时，分别建立10个不同的扩展卡尔曼滤波模型，其中μ_j为针对于不同模型的道路附着系数；这10个模型具有相同的形式，其区别仅在于道路附着系数具体取值的不同；根据三自由度的汽车非线性动力学模型建立扩展卡尔曼滤波的系统状态方程，第j（j=1,2,...10）个模型离散化后的卡尔曼滤波的状态方程的矩阵形式表示为：

X_j(k)＝f_j(X_j(k-1),U_j(k-1),W_j(k-1),γ_j(k-1))               （1）

式（1）中，下标j表示第j个模型（j=1,2,...10），k表示离散化时刻；这10个模型具有相同的系统状态向量，该系统状态向量为X_j＝[x₁ x₂ x₃]′，其中，x₁＝v_x，x₂＝v_y，x₃＝r，v_x、v_y及r分别是汽车的纵向前进速度、侧向速度和横摆角速度，本发明中上角标′表示对矩阵转置；系统外输入向量为U_j＝[u₁ u₂ u₃]′，其中，u₁＝δ，u₂＝F_{j_xf}，u₃＝F_{j_xr}，δ是前轮转向角，F_{j_xf}表示第j个模型（j=1,2,...10）中作用在单个前轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个前轮上的纵向力，F_{j_xr}是第j个模型（j=1,2,...10）中作用在单个后轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个后轮上的纵向力；W_j表示零均值的系统高斯白噪声向量且W_j＝[w₁ w₂ w₃]′，其中w₁、w₂及w₃分别表示三个系统高斯白噪声分量；γ_j表示系统外输入对应的零均值高斯白噪声向量且 ${\mathrm{\γ}}_j={\left[\begin{array}{ccc}{w}_{\mathrm{\δ}}& w_{F_{\mathrm{j\; xf}}}& w_{F_{\mathrm{j\; xr}}}\end{array}\right]}^{\′},$ 其中w_δ表示系统外输入δ对应的零均值高斯白噪声，

和

分别表示外输入F_{j_xf}和F_{j_xr}对应的零均值高斯白噪声，这些白噪声隐含在系统状态函数的系统外输入里面；

非线性的系统状态函数向量为

$f_j(X_j,U_j,W_j,{\mathrm{\γ}}_j)=\left[\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_2}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_3}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\end{array}\right],$

其中，

$\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=v_x(k-1)+\frac{T}{M}[{\mathrm{Mv}}_y(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}\frac{v_y(k-1)+\mathrm{ar}(k-1)}{v_x(k-1)}\mathrm{\δ}(k-1)]+\frac{2T}{M}[F_{j\_\mathrm{xf}}(k-1)+F_{j\_\mathrm{xr}}(k-1)]+w_1\\ f_{j\_2}(X_j(k-1),U_j(k-1),w_j(k-1),{\mathrm{\γ}}_j(k-1))=v_y(k-1)+\frac{T}{M}\{-{\mathrm{Mv}}_x(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}[\mathrm{\δ}(k-1)+\frac{-v_y(k-1)-\mathrm{ar}(k-1)}{v_x(k-1)}]+2C_{\mathrm{\αr}}\frac{\mathrm{br}(k-1)-v_y(k-1)}{v_x(k-1)}\}+\frac{2T}{M}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_2\\ f_{j\_3}(x_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=r(k-1)+\frac{T}{I_z}\left\{2a{C}_{\mathrm{\αf}}\right[\mathrm{\δ}(k-1)-\frac{(v_y(k-1)+\mathrm{ar}(k-1))}{v_x(k-1)}]\\ -2b{C}_{\mathrm{ar}}\frac{[\mathrm{br}(k-1)-v_y(k-1)]}{v_x(k-1)}\}+\frac{2\mathrm{aT}}{I_z}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_3\end{array}$

在上述表达式中，M和I_z分别是车辆的质量和绕过质心垂向轴的转动惯量，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离，C_αf、C_αr分别表示前、后轮胎的侧偏刚度，T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒；W_j对应的系统噪声协方差阵Q_j为：

$Q_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2\end{array}\right],$ 其中

及

分别表示系统高斯白噪声w₁、w₂及w₃对应的方差；γ_j对应的系统外部输入噪声的协方差阵为 ${\mathrm{\Γ}}_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{\mathrm{\δ}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xf}}}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xr}}}}^2\end{array}\right],$ σ_δ ²、

和分别表示w_δ、

和

对应的方差；轮胎纵向力F_{j_xf}和F_{j_xr}根据非线性刷子轮胎模型来确定，轮胎模型中道路附着系数μ_j（j=1,2,...10）取值的不同是10个模型的区别所在；

用s_xq（q=f，r）表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_xf和后轮轴纵向滑移率s_xr，下角标q取f或r，f或r分别表示前或后轮轴，s_xq计算方法为：

s_xq＝(ω_qR-v_xq)/max(ω_qR,v_xq)(q＝f,r)                 （2）

式（2）中，R表示车轮轮胎半径；v_xf和v_xr分别表示前、后轮轴上沿轮胎方向的速度，v_xf和v_xr可统一记为v_xq(q=f,r)；max表示求最大值；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_q(q=f,r)且

$\begin{array}{c}{\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{fR}}+{\mathrm{\ω}}_{\mathrm{fL}})\\ {\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{rR}}+{\mathrm{\ω}}_{\mathrm{rL}})\end{array}----\left(3\right)$

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

v_xq(q=f,r)可按式（4）确定：

$\begin{array}{c}{v}_{\mathrm{xf}}=v_x\cos \mathrm{\δ}+(v_y+\mathrm{ar})\sin \mathrm{\δ}\\ v_{\mathrm{xr}}=v_x\end{array}---\left(4\right)$

进而，轮胎纵向力F_{j_xf}和F_{j_xr}可通过式（5）来确定

$F_{j\_\mathrm{xq}}=C_{\mathrm{xq}}{s}_{\mathrm{xq}}-\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^2}{3{\mathrm{\μ}}_j{F}_{\mathrm{zq}}}+\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^3}{27{\left({\mathrm{\μ}}_j{F}_{\mathrm{zq}}\right)}^2}(q=f,\mathrm{rj}=\mathrm{1,2},...10)---\left(5\right)$

式（5）中，C_xf和C_xr分别表示单个前、后轮胎的纵向刚度，统一记为C_xq(q=f,r)；μ_j（j＝1,2...10）表示轮胎和地面间的道路摩擦系数，所建10个模型的区别仅在于其取值的不同，其中，μ₁＝0.1，μ₂＝0.2...μ₁₀＝1.0；

F_zq(q=f,r)表示分配到单个前或后轮上的垂向载荷且可按下式计算

$F_{\mathrm{zf}}=\frac{\mathrm{Mgb}}{(a+b)},F_{\mathrm{zr}}=\frac{\mathrm{Mga}}{(a+b)}---\left(6\right)$

式（6）中，g表示重力加速度；

车辆纵向前进速度和横摆角速度与两个非转向后轮的速度存在以下关系

$\begin{array}{c}{v}_x=(V_{\mathrm{RL}}+V_{\mathrm{RR}})/2\\ r=(V_{\mathrm{RL}}-V_{\mathrm{RR}})/T_W\end{array}---\left(7\right)$

式（7）中，T_W表示后轮轴上两个后轮间的轮距，V_RL和V_RR分别表示左后轮和右后轮的线速度；

第j（j=1,2,...10）个模型的卡尔曼滤波的观测方程的离散化矩阵形式为：

Z_j(k)=H_j(k)X_j(k)+V_j(k)                   （8）

式（8）中，Z_j为观测向量，H_j为观测阵，V_j表示与W_j互不相关的零均值观测白噪声向量，且 $Z_j\left(k\right)=\left[\begin{array}{c}{v}_{x\_m}\left(k\right)\\ {\mathrm{\ω}}_{z\_m}\left(k\right)\end{array}\right],H_j\left(k\right)=\begin{array}{}\end{array}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],V_j=\left[\begin{array}{c}{n}_{v_x}\\ n_{{\mathrm{\ω}}_z}\end{array}\right],$ 其中v_{x_m}(k)和ω_{z_m}(k)分别为通过轮速传感器测量获得的车辆纵向前进速度和横摆角速度；表示通过轮速传感器测量获得的车辆纵向前进速度的观测噪声且

是均值为0、方差为

的高斯白噪声，

表示通过轮速传感器测量获得的横摆角速度的观测噪声且

是均值为0、方差为

的高斯白噪声；V_j对应的观测噪声方差阵R_j可表示为 $R_j=\left[\begin{array}{cc}{{\mathrm{\σ}}_{v_x}}^2& 0\\ 0& {{\mathrm{\σ}}_{{\mathrm{\ω}}_z}}^2\end{array}\right];$

对于式（8）中的测量值v_{x_m}(k)和ω_{z_m}(k)，它们是利用后轮轴上两个轮速传感器测得的角速度乘以轮胎半径得到V_{RL_m}＝R·ω_rL和V_{RR_m}＝R·ω_rR，V_{RL_m}和V_{RR_m}分别表示V_RL和V_RR含有噪声的测量值，进而利用式（7）获得的，即v_{x_m}和ω_{z_m}分别表示v_x和r的含有噪声的测量值且

2）交互多模型估计方法

对于式（1）描述的系统状态方程和式（8）描述的测量方程，可运用交互多模型滤波理论，建立起滤波递推估计过程。具体估计步骤如下：

①交互估计计算

上述十个扩展卡尔曼系统模型之间的转移概率为p_ij，下标i、j（i=1,2...10,j=1,2,3...10）表示从状态i转移到状态j的概率；

则预测第j(j=1,2...10)个模型的模型概率ρ_j(k,k-1):

${\mathrm{\ρ}}_j=(k,k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}{\mathrm{\ρ}}_i(k-1)$

预测混合概率ρ_i|j(k-1):

ρ_i|j(k-1)＝p_ijρ_i(k-1)/ρ_j(k,k-1)

则交互估计后第j个滤波器的输入为：

$X_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{X}_i(k-1){\mathrm{\ρ}}_{i|j}(k-1)$

$P_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{i|j}(k-1)\{P_i(k-1)+[X_i(k-1)-X_{0j}(k-1)\left]\right[X_i(k-1)-X_{0j}(k-1){]}^{\′}\}$

②每个模型滤波器对于式（1）和式（8）所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，各自进行标准扩展卡尔曼滤波递推，该递推过程包括时间更新和测量更新，第j(j=1,2,...10)个模型的滤波过程如下：

时间更新：

状态一步预测方程X_j(k,k-1)＝f_j(X_0j(k-1),U_j(k-1),0,0)

一步预测误差方差阵：

P_j(k,k-1)＝A_j(k-1)P_0j(k-1)(A_j(k-1))′+B_j(k-1)Γ_j(k-1)(B_j(k-1))′+Q_j(k-1)

其中，A_j、B_j分别是系统状态函数向量f_j对状态向量X_j和外部输入向量U_j求偏导数的雅可比矩阵，即矩阵A_j和B_j的第m行第n列元素A_{j_[m,n]}和B_{j_[m,n]}可分别通过下式求得：

$A_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂x_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

$B_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂u_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

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

$A_{j\_\left[\mathrm{1,1}\right]}=1+T\left[\frac{-2C_{\mathrm{\αf}}(v_y+\mathrm{ar})}{M{v_x}^2}\mathrm{\δ}\right]A_{j\_\left[\mathrm{1,2}\right]}=T[r+\frac{{2C}_{\mathrm{\αf}}}{M_{v_x}}\mathrm{\δ}]$

$A_{j\_\left[\mathrm{1,3}\right]}=T(v_y+\frac{2C_{\mathrm{\αf}}a}{M_{v_x}}\mathrm{\δ})A_{j\_\left[\mathrm{2,1}\right]}=T[-r-2C_{\mathrm{\αr}}\frac{\mathrm{br}-v_y}{M{v_x}^2}+2C_{\mathrm{\αf}}\frac{v_y+\mathrm{ar}}{M{v_x}^2}]$ $A_{j\_\left[\mathrm{2,2}\right]}=1-\frac{2T(C_{\mathrm{\αr}}+C_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}{A}_{j\_\left[\mathrm{2,3}\right]}=T[-v_x+2\frac{(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}]$

$A_{j\_\left[\mathrm{3,1}\right]}=\frac{2T[a{C}_{\mathrm{\αf}}(v_y+\mathrm{ar})+b{C}_{\mathrm{\αr}}(\mathrm{br}-v_y)]}{I_z{v_x}^2}$

$A_{j\_\left[\mathrm{3,2}\right]}=\frac{2T(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{I_z{v}_x}{A}_{j\_\left[\mathrm{3,3}\right]}=1-\frac{2T(a^2{C}_{\mathrm{\αf}}+b^2{C}_{\mathrm{\αr}})}{I_z{v}_x}$

$B_{j\_\left[\mathrm{1,1}\right]}=2T{C}_{\mathrm{\αf}}\frac{(v_y+\mathrm{ar})}{M{v}_x}{B}_{j\_\left[\mathrm{1,2}\right]}=\frac{2T}{M}{B}_{j\_\left[\mathrm{1,3}\right]}=\frac{2T}{M}$

$B_{j\_\left[\mathrm{2,1}\right]}=\frac{2T{F}_{j\_\mathrm{xf}}}{M}+\frac{2T{C}_{\mathrm{\αf}}}{M}{B}_{j\_\left[\mathrm{2,2}\right]}=\frac{2\mathrm{T\δ}}{M}{B}_{j\_\left[\mathrm{2,3}\right]}=0$

$B_{j\_\left[\mathrm{3,1}\right]}=\frac{2\mathrm{Ta}}{I_z}{C}_{\mathrm{\αf}}+\frac{2\mathrm{Ta}}{I_z}{F}_{j\_\mathrm{xf}}{B}_{j\_\left[\mathrm{3,2}\right]}=\frac{2\mathrm{Ta\δ}}{I_z}{B}_{j\_\left[\mathrm{3,3}\right]}=0$

测量更新：

滤波增益矩阵：K_j(k)＝P_j(k,k-1)(H_j(k))′(S_j(k))^-1

S_j(k)＝H_j(k)P_j(k,k-1)(H_j(k))′+R_j(k)

状态估计：X_j(k)＝X_j(k,k-1)+K_j(k)(Z_j(k)-H_j(k)X_j(k,k-1))

估计误差方差阵：P_j(k)＝P_j(k,k-1)-K_j(k)S_j(k)(K_j(k))′

③模型概率更新

在每个模型完成上一步的更新之后，利用最大似然函数Λ_j(k)计算新的模型概率ρ_j(k)，最大似然函数计算如下：

${\mathrm{\Λ}}_j\left(k\right)=\frac{\exp \{-\frac{1}{2}{(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))}^{\′}{\left(S_j\left(k\right)\right)}^{-1}(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))\}}{|2\mathrm{\π}{S}_j\left(k\right){|}^{-\frac{1}{2}}}$

因此，模型j在k时刻的模型概率由贝叶斯定理给出：

${\mathrm{\ρ}}_j\left(k\right)=\frac{{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_j(k,k-1)}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_i(k,k-1)}$

④估计组合

在计算出各模型为正确的后验概率之后，对所有滤波器的状态估计进行概率加权并求和，权系数为模型正确的后验概率，得到最终的状态估计为：

$X\left(k\right)=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{X}_j\left(k\right){\mathrm{\ρ}}_j\left(k\right),$ 其中， $X\left(k\right)={\left[\begin{array}{ccc}\stackrel{\‾}{v_x}& \stackrel{\‾}{v_y}& \stackrel{\‾}{r}\end{array}\right]}^{\′},$ 各状态变量的上标“-”表示各状态量的最终滤波估计值，即X(k)内各状态变量依次分别表示估计组合后的纵向车速、侧向车速和横摆角速度；

同时，由于各模型的区别在于各模型所设定的道路附着系数的具体取值不同，即各模型的μ_j的取值不同，因此，对各模型所设定的附着系数进行概率加权即可得出最终滤波估计出的当前的道路附着系数μ：

$\mathrm{\μ}=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\μ}}_j{\mathrm{\ρ}}_j\left(k\right)---\left(9\right).$

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

1.本发明的方法是针对不同附着系数条件，在非线性整车动力学模型和多个轮胎纵向力模型基础上提出的，在不同道路附着系数条件下仍可以获得准确的车辆纵向和侧向车速信息，同时实时给出相应的道路附着系数估计，可用于汽车主动安全控制对车辆纵向、横向车速以及道路附着系数的测量与估计需要，具有精度高、成本低、实时性好、全面自适应等优点。

2.本发明提出的基于交互多模型的车速与道路附着系数的联合估计方法对于道路附着系数的突变具有良好的适应性，响应时间短，能够满足汽车纵向主动安全控制的要求。

3.本发明提出的基于交互多模型的车速与道路附着系数的联合估计方法采用3自由度非线性模型对车辆的纵向前进速度并无定常的限定，故即可适应一般机动环境也可适应较高机动环境下车辆运行状态的准确估计，达到全面适应的效果。

附图说明

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

图2.车辆动力学模型

图3.垂直载荷一定时，刷子轮胎模型纵向力与滑移率的关系

图4.单附着系数路面仿真设定的方向盘转角（度）随时间变化图

图5.单附着系数路面仿真设定的纵向速度（米/秒）随时间变化图

图6.单附着系数路面本发明方法对道路附着系数估计结果（图中估计结果用

灰虚线代表，Carsim输出真值用黑实线代表）

图7.单附着系数路面本发明方法对纵向车速估计误差

图8.单附着系数路面普通扩展卡尔曼滤波方法对纵向车速估计误差

图9.附着系数突变路面仿真设定的方向盘转角（度）随时间变化图

图10.附着系数突变路面仿真设定的纵向速度（米/秒）随时间变化图

图11.附着系数突变路面本发明方法对道路附着系数估计结果（图中估计结

果用灰虚线代表，Carsim输出真值用黑实线代表）

具体实施方式

实施实例1

随着社会经济的发展，道路交通安全问题日益突出，并已成为全球性难题。全世界每年因交通事故都会造成大量的人员伤亡和财产损失，世界各国都在努力降低交通事故的发生。近年来，汽车主动安全技术得到了迅速的发展。汽车主动安全技术能够防患于未然，主动避免事故的发生，已成为现代汽车最主要的发展方向之一。目前常见的主动安全技术主要包括防抱死制动系统（ABS），车辆电子稳定程序（ESP），牵引力控制系统（TCS），电控驱动防滑系统（ASR），四轮转向稳定控制系统（4WS）等。这些系统通常涉及汽车的纵向前进速度、侧向速度、以及质心侧偏角等运行状态的测量或估计，而这些运行状态的测量可用于后续的汽车主动安全控制，因此其精度直接关系汽车的行驶安全性与稳定性，追尾碰撞预警/避免系统（CW/CA）、制动防抱死系统（ABS）等纵向主动安全系统依赖于纵向车速的准确估计，而以电子稳定程序（ESP）、四轮转向控制（4WS）为代表的侧向安全系统则依赖于对质心侧偏角的准确估计，而质心侧偏角可根据纵向车速和侧向车速计算得知。所以说，纵向和侧向车速信息获取的准确性和可靠性，直接决定着这些主动安全系统的有效性。同时，这些主动安全系统效果的优劣很大程度上取决于能否“道路自适应”，即如果能够实时估计出道路附着系数，系统就可以根据当前路况调节控制策略，提高车辆安全。因此，纵向、侧向车速与道路附着系数作为汽车主动安全的关键参数，对其进行实时、准确的测量或估计，是上述汽车主动安全系统发挥作用的重要前提与基础。

但是由于成本、技术等方面的原因，这些车速和道路附着系数信息通常难以直接测量，而是使用量产车上已装备的传感器，通过建立车辆运动或动力学模型，利用车辆状态估计技术来获得。

在车速估计方面，主要有基于运动学模型和基于动力学模型的估计方法，其中，基于运动学模型的估计方法虽然鲁棒性较好，不受模型参数影响，但对于传感器精度要求较高，车载低成本传感器的精度往往难以满足其需求，给其应用带来了很大的限制。而基于动力学模型的估计方法则对传感器精度要求相对较低，适宜于与车载传感器结合使用，因此成为了目前较为普遍的方法。在三自由度整车动力学模型或两自由度线性整车模型的基础上，结合线性或各种非线性轮胎模型，卡尔曼滤波算法、滑模观测器,模糊观测器,神经网络观测器及各种非线性观测器被广泛应用于车速估计中，取得了较好的效果。但基于动力学模型的方法对于模型本身的精度要求很高，要求模型尽可能准确的反应车辆的动力学特性，因而对模型参数准确性要求较高，同时对参数的变化较为敏感，道路附着系数作为重要的轮胎力学参数，在这些方法中往往被假设为已知且定常，从而根据事先设定的附着系数经验值，建立单一的轮胎模型计算轮胎力以估计车速，事实上，道路附着系数往往偏离事先预设的经验值，且在车辆运行过程中会变化而非定常，从而使模型不能很好的适应于当前道路情况，或导致车速估计结果无法及时响应道路条件的改变，严重影响其准确性。且目前已提出的动力学模型对整车或轮胎做了较多线性化假定的动力学模型，这些模型在车辆较平稳运行时能获得较好的估计效果和精度，但在较高机动运行状况下由于难于反映车辆的实际非线性动力学行为导致估计精度较低。

目前，在汽车主动安全领域，道路附着系数主要分为直接测量和间接估计两类，直接测量方法是利用光、声、图像、雷达等传感器直接检测路面，测量一些对路面附着系数影响较大的因素，并根据以往经验预测当前道路附着系数的大小，但这些方法都需要额外加装传感器，且传感器成本都较高，难以实现大规模的商业应用，其次需要进行大量的测试训练，识别精度很大程度上依赖于经验，难以准确估算没有测试和训练过的路况的附着系数。间接估计方法是通过对汽车的运行过程进行运动学或动力学建模，结合轮胎模型，将有关低成本的车载传感器（如轮速传感器、陀螺仪、加速度计以及GPS等）信息作为观测信息，进而利用适当的滤波估计算法实现对道路附着系数的估计。已有的间接方法包括基于车辆侧向动力学和基于纵向动力学的研究两种，但这些估计方法或有其特定的适用范围（例如仅适宜于滑移率较低的线性区域），或需要一定的经验性，或存在难以测量的参数，或在实际使用中的效果有待进一步确认，同时，道路附着系数的准确估计也依赖于对于纵向、横向车速等车辆运行状态的准确估计，这些都限制了这些估计方法在汽车主动安全系统上的应用。

总而言之，由于道路附着系数和车辆纵向、横向车速在车辆运行过程中相互作用，相互影响，传统估计方法中对于道路附着系数定常的假设无法满足不同道路附着系数条件下对车速的准确估计的要求，也无法实时估计出道路附着系数值。

为实现在不同道路附着系数条件下对车辆纵向、侧向车速的准确估计，以满足汽车主动安全系统的需求，并实时估计出道路附着系数值，本发明提出了一种基于交互多模型（Interacting Multiple Model，IMM）的车速与道路附着系数的联合估计方法。本发明提出的估计方法可在不同道路附着系数条件下实现对车辆纵向车速和侧向车速的准确估计，同时实时给出相应的道路附着系数估计，具有精度高、成本低、实时性好、全面自适应等特点，本发明的具体思路如下：

交互多模型算法具有自适应的特点，通过建立不同的多个模型滤波器，各模型滤波器通过估计状态的组合实现交互，模型之间基于马尔可夫链进行切换，能够有效地对各个模型的概率进行调整。本发明的交互多模型算法中通过扩展卡尔曼滤波(Extended Kalman Filter，EKF)方法建立模型滤波器，卡尔曼滤波器是以最小均方差为准则的最优状态估计滤波器，它不需要储存过去的测量值，只根据当前的观测值和前一时刻的估计值，利用计算机进行递推计算，便可实现对实时信号的估计。递归最小二乘和卡尔曼滤波都具有数据存储量小、算法简便的特点。

为适应不同道路附着系数环境和较高机动环境下汽车主动安全控制对车辆纵向、侧向车速信号以及道路附着系数的测量与估计要求，首先对汽车以及轮胎进行适当的动力学建模，即建立卡尔曼滤波过程的系统状态方程。针对本发明的应用领域，本发明对于行驶在通常道路交通环境上的前轮转向的四轮车辆(目前应有最广的情况，典型例子如前轮转向的轿车)，可做如下的合理假定：

1)忽略汽车的俯仰、侧倾和上下弹跳运动。

2)忽略汽车悬架对轮胎轴上的影响。

3)忽略侧倾运动，可认为汽车前轴上左右两个轮胎的转向角、侧偏角、纵向力及侧向力相同；类似地，可假定汽车后轴上左右两个轮胎的侧偏角、纵向力及侧向力相同。

根据上述应用要求和假定，本发明针对目前应用较多的前轮转向四轮汽车，采用附图2所示的车辆动力学模型(经等效简化后相当于前、后车轮被分别集中在汽车前、后轴中点而构成的一假想Bicycle模型，如图2右侧所示)。该模型有3个自由度，分别是纵向运动、侧向运动以及横摆转动。图2中定义了车辆载体坐标系，其原点o位于质心处，ox轴沿车辆的纵向轴并与车辆前进方向一致，oz轴垂直于车辆运行平面并指向地面(即向下，绕oz轴的横摆角速度r的正方向定义如图示)，而oy轴按右手螺旋规则可确定。纵向前进速度v_x、侧向速度v_y和横摆角速度r都是指车辆质心的。根据牛顿力学，车辆的动力学模型可描述为

纵向： $\begin{array}{cc}{\stackrel{\·}{v}}_x=v_y\·r+a_x& \left(1\right)\\ a_x=\frac{2}{M}[F_{\mathrm{xf}}\cos \left(\mathrm{\δ}\right)-F_{\mathrm{yf}}\sin \left(\mathrm{\δ}\right)+F_{\mathrm{xr}}& \left(2\right)\end{array}$

横向： $\begin{array}{cc}{\stackrel{\·}{v}}_y=-v_x\·r+a_y& \left(3\right)\\ a_y=\frac{2}{M}[F_{\mathrm{xf}}\sin \left(\mathrm{\δ}\right)+F_{\mathrm{yf}}\cos \left(\mathrm{\δ}\right)+F_{\mathrm{yr}}& \left(4\right)\end{array}$

横摆：

$\stackrel{\·}{r}=\frac{2}{I_z}[a{F}_{\mathrm{xf}}\sin \left(\mathrm{\δ}\right)+a\·F_{\mathrm{yf}}\cos \left(\mathrm{\δ}\right)-b{F}_{\mathrm{yr}}]---\left(5\right)$

式中，v_x、v_y及r分别是汽车的纵向前进速度、侧向速度和横摆角速度，本发明中，上标志“·”表示微分，如

表示对r的微分；a_x与a_y分别是汽车纵向和侧向加速度；M和I_z分别是车辆的质量和绕oz轴的转动惯量；a、b分别是汽车前、后轮轮轴中心到质心的距离；δ是前轮转向角；F_xf和F_xr是作用在单个前轮和后轮上的纵向力；F_yf和F_yr是作用在单个前轮和后轮上的侧向力。

对于在一般道路行驶的车辆，通常可将作用在各轮上的侧向力表示为：

F_yf＝C_αfα_f，F_yr＝C_αrα_r                      （6）

式（6）中，C_αf、C_αr分别是前、后轮胎的侧偏刚度；α_f、α_r分别是前、后轮胎的侧偏角且可表示为

${\mathrm{\α}}_f=\mathrm{\δ}-\frac{v_y+\mathrm{ar}}{v_x},{\mathrm{\α}}_r=\frac{\mathrm{br}-v_y}{v_x}---\left(7\right)$

将式（6）、（7）代入式（1）-（5），并考虑到δ通常是小角度，即cos（δ）≈1，sin（δ）≈δ；且忽略二阶及以上的高阶微量，经整理后可得：

$\left\{\begin{array}{c}{\stackrel{\·}{v}}_x=\frac{1}{M}[M{v}_{y}r+2\frac{v_y+\mathrm{ar}}{v_x}{C}_{\mathrm{\αf}}\mathrm{\δ}]+\frac{2}{M}(F_{\mathrm{xf}}+F_{\mathrm{xr}})\\ {\stackrel{\·}{v}}_y=\frac{1}{M}[-M{v}_{x}r+2(\mathrm{\δ}-\frac{(v_y+\mathrm{ar})}{v_x})C_{\mathrm{\αf}}+2C_{\mathrm{\αr}}\frac{\mathrm{br}-v_y}{v_x}]+\frac{2}{M}{F}_{\mathrm{xf}}\mathrm{\δ}\\ \stackrel{\·}{r}=\frac{1}{I_z}[2a(\mathrm{\δ}-\frac{(v_y+\mathrm{ar})}{v_x})C_{\mathrm{\αf}}-2b{C}_{\mathrm{\αr}}\frac{(\mathrm{br}-v_y)}{v_x}]+\frac{2a}{I_z}{F}_{\mathrm{xf}}\mathrm{\δ}\end{array}\right.---\left(8\right)$

对于式（8）中的前轮转向角δ可通过方向盘转角传感器测得的方向盘转角除以从方向盘到前轮的转向传动比来确定。而对于式（8）中的轮胎纵向力，采用轮胎模型来确定。魔术公式轮胎模型是公认的拟合精度最高的经验轮胎模型，但它是由三角函数组合而成的复杂的非线性函数，且模型中未知因子较多，计算量较大，不适于实时使用。因此，本发明中纵向力的估计确定采用便于实时计算的非线性刷子模型[可参考文献:Pacejka H B.Analysis of tire properties.In:Clark SK(ed.).Mechanics of Pneumatic Tires,new edition.Washington DC:DOT HS805952,NHTSA,1981:721～870]。为此，引入车辆纵向滑移率s_xq（q=f，r）表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_xf和后轮轴纵向滑移率s_xr，本发明中下角标q取f或r，f或r分别表示前或后轮轴，s_xq计算方法为：

s_xq＝(ω_qR-v_xq)/max(ω_qR,v_xq)(q＝f,r)            （9）

式（9）中，R表示车轮轮胎半径；v_xf和v_xr分别表示前、后轮轴上沿轮胎方向的速度，v_xf和v_xr可统一记为v_xq(q=f,r)；max表示求最大值运算；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_q(q=f,r)且

$\begin{array}{c}{\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{fR}}+{\mathrm{\ω}}_{\mathrm{fL}})\\ {\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{rR}}+{\mathrm{\ω}}_{\mathrm{rL}})\end{array}---\left(10\right)$

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

v_xq(q=f,r)可按式（11）确定：

$\begin{array}{c}{v}_{\mathrm{xf}}=v_x\cos \mathrm{\δ}+(v_y+\mathrm{ar})\sin \mathrm{\δ}\\ v_{\mathrm{xr}}=v_x\end{array}---\left(11\right)$

由于本发明采取交互多模型方法，对于道路附着系数分别为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0建立10个不同的扩展卡尔曼滤波模型，用μ_j（j=1,2...10）表示第j（j=1,2,...10）个模型中轮胎和地面间的道路摩擦系数，其中，μ_j＝10×j（j=1,2,...10），即μ₁＝0.1，μ₂＝0.2...μ₁₀＝1.0，对应于每个模型，作用在单个前轮和后轮的纵向力分别用F_{j_xf}和F_{j_xr}（j=1,2,...10）表示，F_{j_xf}和F_{j_xr}可统一记为F_{j_xq}（j=1,2,...10），下角标q取f或r。则轮胎纵向力可通过式（12）来确定：

$F_{j\_\mathrm{xq}}=C_{\mathrm{xq}}{s}_{\mathrm{xq}}-\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^2}{3{\mathrm{\μ}}_j{F}_{\mathrm{zq}}}+\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^3}{27{\left({\mathrm{\μ}}_j{F}_{\mathrm{zq}}\right)}^2}(q=f,\mathrm{rj}=\mathrm{1,2},...10)---\left(12\right)$

式（12）中，C_xf和C_xr分别表示单个前、后轮胎的纵向刚度，统一记为C_xq(q=f,r)；本发明所建立的多模型即是针对道路附着系数分别取不同值时所对应的模型，即本发明中所建立的多模型具有相同的形式，其区别仅仅在于道路附着系数具体取值的不同（取值为0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0十个）即μ_j＝10×j（j=1,2,...10），μ_j为针对于不同模型的道路附着系数，附图3表示了垂直载荷一定时，道路附着系数取值分别等于0.1,0.2,…1时，纵向力与滑移率之间的关系。而在传统扩展卡尔曼滤波方法中，道路附着系数取值是一成不变的，为预先预设的经验值，当预设值与实际道路附着系数相差较大或道路附着系数发生突变时，往往导致车速估计结果误差很大。

而F_zq(q=f,r)表示分配到单个前或后轮胎上的垂向载荷且可按下式计算

$F_{\mathrm{zf}}=\frac{\mathrm{Mgb}}{(a+b)},F_{\mathrm{zr}}=\frac{\mathrm{Mga}}{(a+b)}---\left(13\right)$

式（13）中，g表示重力加速度。

车辆纵向前进速度和横摆角速度与两个非转向后轮的速度存在以下关系

$\begin{array}{c}{v}_x=(V_{\mathrm{RL}}+V_{\mathrm{RR}})/2\\ r=(V_{\mathrm{RL}}-V_{\mathrm{RR}})/T_W\end{array}---\left(7\right)$

式（14）中，T_W表示后轮轴上两个后轮间的轮距，V_RL和V_RR分别表示左后轮和右后轮的线速度。

对于式（8）描述的模型，它是一个具有3自由度的非线性车辆动力学模型，不同于经常所采用的2自由度线性车辆模型。在经常采用的2自由度线性车辆模型中，车辆的纵向前进速度被认为是定常的，车辆模型仅是关于侧向速度和横摆角速度的线性微分方程。因此，2自由度线性车辆模型一般只适合前向速度不变或变化缓慢的运行情况（机动性较低），而对于较高机动运行情况（即需要频繁转向以及加减速的情形），该模型存在较大的建模误差。而本发明所采用的3自由度非线性模型对车辆的纵向前进速度并无定常的限定，故即可适应一般机动环境也可适应较高机动环境下车辆运行状态的准确估计。因此，本发明将根据式（8）建立卡尔曼滤波的系统状态方程。

本发明采用的交互多模型算法中将设置多个滤波模型，这多个模型是针对道路附着系数分别取值为0.1、0.2....1.0，即针对不同模型的道路附着系数μ_j（j=1,2,...10）分别取值为0.1、0.2....1.0时，所分别建立的10个不同的扩展卡尔曼滤波模型，因此，所建立的卡尔曼滤波的状态方程也应有10个。而这10个模型具有相同的形式，其区别仅在于道路附着系数具体取值的不同。系统模型之间的转移概率为p_ij，下标i、j（i=1,2,...10,j=1,2,...10）表示从状态i转移到状态j的概率；

应注意的是，在实际的卡尔曼滤波递推过程中，需采用离散化的卡尔曼滤波模型。为此，对式（8）的微分方程组进行离散化处理，这十个模型具有相同的形式，其区别仅在于道路附着系数具体取值的不同，第j个EKF模型建立离散化后的系统方程和观测方程为：

X_j(k)＝f_j(X_j(k-1),U_j(k-1),W_j(k-1),γ_j(k-1))            （15）

式（15）中，下标j表示第j个模型（j=1,2,...10），k表示离散化时刻；这10个模型具有相同的系统状态向量，该系统状态向量为X_j＝[x₁ x₂ x₃]′，其中，x₁＝v_x，x₂＝v_y，x₃＝r，v_x、v_y及r分别是汽车的纵向前进速度、侧向速度和横摆角速度，本发明中上角标′表示对矩阵转置；系统外输入向量为U_j＝[u₁ u₂ u₃]′，其中，u₁＝δ，u₂＝F_{j_xf}，u₃＝F_{j_xr}，δ是前轮转向角，F_{j_xf}表示第j个模型（j=1,2,...10）中作用在单个前轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个前轮上的纵向力，F_{j_xr}是第j个模型（j=1,2,...10）中作用在单个后轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个后轮上的纵向力；W_j表示零均值的系统高斯白噪声向量且W_j＝[w₁ w₂ w₃]′，其中w₁、w₂及w₃分别表示三个系统高斯白噪声分量；γ_j表示系统外输入对应的零均值高斯白噪声向量且 ${\mathrm{\γ}}_j={\left[\begin{array}{ccc}{w}_{\mathrm{\δ}}& w_{F_{\mathrm{j\; xf}}}& w_{F_{\mathrm{j\; xr}}}\end{array}\right]}^{\′},$ 其中w_δ表示系统外输入δ对应的零均值高斯白噪声，

和

分别表示外输入F_{j_xf}和F_{j_xr}对应的零均值高斯白噪声，这些白噪声隐含在系统状态函数的系统外输入里面；

非线性的系统状态函数向量为

$f_j(X_j,U_j,W_j,{\mathrm{\γ}}_j)=\left[\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_2}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_3}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\end{array}\right],$

其中，

$\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=v_x(k-1)+\frac{T}{M}[{\mathrm{Mv}}_y(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}\frac{v_y(k-1)+\mathrm{ar}(k-1)}{v_x(k-1)}\mathrm{\δ}(k-1)]+\frac{2T}{M}[F_{j\_\mathrm{xf}}(k-1)+F_{j\_\mathrm{xr}}(k-1)]+w_1\\ f_{j\_2}(X_j(k-1),U_j(k-1),w_j(k-1),{\mathrm{\γ}}_j(k-1))=v_y(k-1)+\frac{T}{M}\{-{\mathrm{Mv}}_x(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}[\mathrm{\δ}(k-1)+\frac{-v_y(k-1)-\mathrm{ar}(k-1)}{v_x(k-1)}]+2C_{\mathrm{\αr}}\frac{\mathrm{br}(k-1)-v_y(k-1)}{v_x(k-1)}\}+\frac{2T}{M}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_2\\ f_{j\_3}(x_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=r(k-1)+\frac{T}{I_z}\left\{2a{C}_{\mathrm{\αf}}\right[\mathrm{\δ}(k-1)-\frac{(v_y(k-1)+\mathrm{ar}(k-1))}{v_x(k-1)}]\\ -2b{C}_{\mathrm{ar}}\frac{[\mathrm{br}(k-1)-v_y(k-1)]}{v_x(k-1)}\}+\frac{2\mathrm{aT}}{I_z}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_3\end{array}$

在上述表达式中，M和I_z分别是车辆的质量和绕过质心垂向轴的转动惯量，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离，C_αf、C_αr分别表示前、后轮胎的侧偏刚度，T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒；W_j对应的系统噪声协方差阵Q_j为：

$Q_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2\end{array}\right],$ 其中及分别表示系统高斯白噪声w₁、w₂及w₃对应的方差；γ_j对应的系统外部输入噪声的协方差阵为 ${\mathrm{\Γ}}_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{\mathrm{\δ}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xf}}}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xr}}}}^2\end{array}\right],$

和

分别表示w_δ、

和

对应的方差；轮胎纵向力F_{j_xf}和F_{j_xr}根据非线性刷子轮胎模型来确定，轮胎模型中道路附着系数μ_j（j=1,2,...10）取值的不同是10个模型的区别所在；

建立车辆运行状态估计的卡尔曼滤波模型的系统状态方程后，下面讨论如何建立其观测方程。从运动学角度，图2所示的车辆运动实际上是一个平面复合运动（纵向运动、侧向运动和横摆转动的复合），故根据平面复合运动关系，可得

$\begin{array}{c}{V}_{\mathrm{RL}}=v_x+\frac{T_W}{2}r\\ V_{\mathrm{RR}}=v_x-\frac{T_W}{2}r\end{array}---\left(16\right)$

式中，V_RL和V_RR分别代表左后轮和右后轮（即两个非转向轮）的车轮线速度，T_W是后轮轴上两个后轮间的轮距。

对式（16）重新整理，可以得到

$\begin{array}{c}{v}_x=(V_{\mathrm{RL}}+V_{\mathrm{RR}})/2\\ r=(V_{\mathrm{RL}}-V_{\mathrm{RR}})/T_W\end{array}---\left(17\right)$

需要指出的是，左后轮和右后轮的车轮线速度可通过安装在后轮轴上的两个轮速传感器获得，即利用后轮轴上两个轮速传感器测得的角速度乘以轮胎半径得到。考虑到轮速传感器的测量噪声，V_{RL_m}＝R·ω_rL与V_{RR_m}＝R·ω_rR，其中V_{RL_m}和V_{RR_m}分别表示V_RL和V_RR含有噪声的测量值。另外，V_{RL_m}和V_{RR_m}还可分别表示为 $V_{\mathrm{RL}\_m}=V_{\mathrm{RL}}+n_{V_{\mathrm{RL}}},V_{\mathrm{RR}\_m}=V_{\mathrm{RR}}+n_{V_{\mathrm{RR}}},$ 其中

和

分别表示左后轮和右后轮的车轮线速度的加性测量噪声（均可建模为均值为0的高斯白噪声）。

在本发明中，将纵向前进速度和横摆角速度作为卡尔曼滤波模型的观测量。由于纵向前进速度和横摆角速度同时又是上述建立的卡尔曼滤波模型的两个状态，故不难建立滤波系统的观测方程，第j（j=1,2,...10）个模型的卡尔曼滤波的观测方程的离散化矩阵形式为：

Z_j(k)=H_j(k)X_j(k)+V_j(k)

      （18）

式（18）中，Z_j为观测向量，H_j为观测阵，V_j表示与W_j互不相关的零均值观测白噪声向量，且 $Z_j\left(k\right)=\left[\begin{array}{c}{v}_{x\_m}\left(k\right)\\ {\mathrm{\ω}}_{z\_m}\left(k\right)\end{array}\right],H_j\left(k\right)=\begin{array}{}\end{array}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],V_j=\left[\begin{array}{c}{n}_{v_x}\\ n_{{\mathrm{\ω}}_z}\end{array}\right],$ 其中v_{x_m}(k)和ω_{z_m}(k)分别为通过轮速传感器测量获得的车辆纵向前进速度和横摆角速度；

表示通过轮速传感器测量获得的车辆纵向前进速度的观测噪声且

是均值为0、方差为

的高斯白噪声，

表示通过轮速传感器测量获得的横摆角速度的观测噪声且

是均值为0、方差为

的高斯白噪声；V_j对应的观测噪声方差阵R_j可表示为 $R_j=\left[\begin{array}{cc}{{\mathrm{\σ}}_{v_x}}^2& 0\\ 0& {{\mathrm{\σ}}_{{\mathrm{\ω}}_z}}^2\end{array}\right];$

对于式（18）中的测量值v_{x_m}(k)和ω_{z_m}(k)，它们是利用后轮轴上两个轮速传感器测得的角速度乘以轮胎半径得到V_{RL_m}＝R·ω_rL和V_{RR_m}＝R·ω_rR，V_{RL_m}和V_{RR_m}分别表示V_RL和V_RR含有噪声的测量值，进而利用式（17）获得的，即v_{x_m}和ω_{z_m}分别表示v_x和r的含有噪声的测量值且

对于式（15）描述的系统状态方程和式（18）描述的测量方程，可运用交互多模型滤波理论，建立起滤波递推估计过程。具体估计步骤如下：

（1）交互估计计算

上述十个扩展卡尔曼系统模型之间的转移概率为p_ij，下标i、j（i=1,2...10,j=1,2,3...10）表示从状态i转移到状态j的概率；

则预测第j(j=1,2...10)个模型的模型概率ρ_j(k,k-1):

${\mathrm{\ρ}}_j=(k,k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}{\mathrm{\ρ}}_i(k-1)$

预测混合概率ρ_i|j(k-1):

ρ_i|j(k-1)＝p_ijρ_i(k-1)/ρ_j(k,k-1)

则交互估计后第j个滤波器的输入为：

$X_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{X}_i(k-1){\mathrm{\ρ}}_{i|j}(k-1)$

$P_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{i|j}(k-1)\{P_i(k-1)+[X_i(k-1)-X_{0j}(k-1)\left]\right[X_i(k-1)-X_{0j}(k-1){]}^{\′}\}$

（2）每个模型滤波器对于式（15）和式（18）所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，各自进行标准扩展卡尔曼滤波递推，但注意到式（15）所示的状态方程为非线性方程，在应用卡尔曼滤波计算时，需先进行线性化处理，将系统方程在X(k,k-1)附近按泰勒级数展开，保留一阶微量、忽略高阶微量后再进行滤波递推计算，即需按照扩展卡尔曼滤波过程进行滤波递推。该递推过程包括时间更新和测量更新，第j(j=1,2,3...10)个模型的滤波过程如下：

时间更新：

状态一步预测方程X_j(k,k-1)＝f_j(X_0j(k-1),U_j(k-1),0,0)

一步预测误差方差阵：

P_j(k,k-1)＝A_j(k-1)P_0j(k-1)(A_j(k-1))′+B_j(k-1)Γ_j(k-1)(B_j(k-1))′+Q_j(k-1)

其中，A_j、B_j分别是系统状态函数向量f_j对状态向量X_j和外部输入向量U_j求偏导数的雅可比矩阵，即矩阵A_j和B_j的第m行第n列元素A_{j_[m,n]}和B_{j_[m,n]}可分别通过下式求得：

$A_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂x_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

$B_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂u_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

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

$A_{j\_\left[\mathrm{1,1}\right]}=1+T\left[\frac{-2C_{\mathrm{\αf}}(v_y+\mathrm{ar})}{M{v_x}^2}\mathrm{\δ}\right]A_{j\_\left[\mathrm{1,2}\right]}=T[r+\frac{{2C}_{\mathrm{\αf}}}{M_{v_x}}\mathrm{\δ}]$

$A_{j\_\left[\mathrm{1,3}\right]}=T(v_y+\frac{2C_{\mathrm{\αf}}a}{M_{v_x}}\mathrm{\δ})A_{j\_\left[\mathrm{2,1}\right]}=T[-r-2C_{\mathrm{\αr}}\frac{\mathrm{br}-v_y}{M{v_x}^2}+2C_{\mathrm{\αf}}\frac{v_y+\mathrm{ar}}{M{v_x}^2}]$ $A_{j\_\left[\mathrm{2,2}\right]}=1-\frac{2T(C_{\mathrm{\αr}}+C_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}{A}_{j\_\left[\mathrm{2,3}\right]}=T[-v_x+2\frac{(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}]$

$A_{j\_\left[\mathrm{3,1}\right]}=\frac{2T[a{C}_{\mathrm{\αf}}(v_y+\mathrm{ar})+b{C}_{\mathrm{\αr}}(\mathrm{br}-v_y)]}{I_z{v_x}^2}$

$A_{j\_\left[\mathrm{3,2}\right]}=\frac{2T(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{I_z{v}_x}{A}_{j\_\left[\mathrm{3,3}\right]}=1-\frac{2T(a^2{C}_{\mathrm{\αf}}+b^2{C}_{\mathrm{\αr}})}{I_z{v}_x}$

$B_{j\_\left[\mathrm{1,1}\right]}=2T{C}_{\mathrm{\αf}}\frac{(v_y+\mathrm{ar})}{M{v}_x}{B}_{j\_\left[\mathrm{1,2}\right]}=\frac{2T}{M}{B}_{j\_\left[\mathrm{1,3}\right]}=\frac{2T}{M}$

$B_{j\_\left[\mathrm{2,1}\right]}=\frac{2T{F}_{j\_\mathrm{xf}}}{M}+\frac{2T{C}_{\mathrm{\αf}}}{M}{B}_{j\_\left[\mathrm{2,2}\right]}=\frac{2\mathrm{T\δ}}{M}{B}_{j\_\left[\mathrm{2,3}\right]}=0$

$B_{j\_\left[\mathrm{3,1}\right]}=\frac{2\mathrm{Ta}}{I_z}{C}_{\mathrm{\αf}}+\frac{2\mathrm{Ta}}{I_z}{F}_{j\_\mathrm{xf}}{B}_{j\_\left[\mathrm{3,2}\right]}=\frac{2\mathrm{Ta\δ}}{I_z}{B}_{j\_\left[\mathrm{3,3}\right]}=0$

测量更新：

滤波增益矩阵：K_j(k)＝P_j(k,k-1)(H_j(k))′(S_j(k))^-1

S_j(k)＝H_j(k)P_j(k,k-1)(H_j(k))′+R_j(k)

状态估计：X_j(k)＝X_j(k,k-1)+K_j(k)(Z_j(k)-H_j(k)X_j(k,k-1))

估计误差方差阵：P_j(k)＝P_j(k,k-1)-K_j(k)S_j(k)(K_j(k))′

（3）模型概率更新

在每个模型完成上一步的更新之后，利用最大似然函数Λ_j(k)计算新的模型概率ρ_j(k)，最大似然函数计算如下：

${\mathrm{\Λ}}_j\left(k\right)=\frac{\exp \{-\frac{1}{2}{(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))}^{\′}{\left(S_j\left(k\right)\right)}^{-1}(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))\}}{|2\mathrm{\π}{S}_j\left(k\right){|}^{-\frac{1}{2}}}$

因此，模型j在k时刻的模型概率由贝叶斯定理给出：

${\mathrm{\ρ}}_j\left(k\right)=\frac{{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_j(k,k-1)}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_i(k,k-1)}$

（4）估计组合

在计算出各模型为正确的后验概率之后，对所有滤波器的状态估计进行概率加权并求和，权系数为模型正确的后验概率，得到最终的状态估计为：

$X\left(k\right)=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{X}_j\left(k\right){\mathrm{\ρ}}_j\left(k\right),$ 其中， $X\left(k\right)={\left[\begin{array}{ccc}\stackrel{\‾}{v_x}& \stackrel{\‾}{v_y}& \stackrel{\‾}{r}\end{array}\right]}^{\′},$ 各状态变量的上标“-”表示各状态量的最终滤波估计值，即X(k)内各状态变量依次分别表示估计组合后的纵向车速、侧向车速和横摆角速度；

同时，由于各模型的区别在于各模型所设定的道路附着系数的具体取值不同，即各模型的μ_j的取值不同，因此，对各模型所设定的附着系数进行概率加权即可得出最终滤波估计出的当前的道路附着系数μ：

$\mathrm{\μ}=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\μ}}_j{\mathrm{\ρ}}_j\left(k\right)---\left(19\right)$

实施实例2

为检验本发明提出的车速与道路附着系数的联合估计方法的实际效果，在专业的汽车动力学仿真软件CarSim上进行了仿真验证实验。

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

仿真所用车辆是一个前轮转向的四轮车，主要参数如下：M=960（千克）、Iz=1382（千克·米·米）、a=0.948（米）、b=1.422（米）、C_αf=C_αr=25927（牛顿/弧度）、T_w=1.390（米）。设定四个车轮的线速度（通过轮速传感器测得的角速度乘以轮胎半径得到）的测量噪声均为均值是0、标准差是0.04（米/秒）的高斯白噪声，方向盘转角传感器的测量噪声为均值是0、标准差是0.0873（弧度）的高斯白噪声。卡尔曼滤波的系统零均值高斯白噪声的标准差分别为

及

卡尔曼滤波的外输入的零均值高斯白噪声的标准差为σ_δ＝0.00873（弧度），

及

卡尔曼滤波的两个观测量的零均值高斯白噪声的标准差分别为

（米/秒）及

为检验本发明所提出估计方法对于不同路面状况的适应性，分别针对单附着系数路面与附着系数跃变的路面对本发明提出算法进行验证，并与传统的扩展卡尔曼滤波（EKF）算法所估计结果进行比较，传统的扩展卡尔曼滤波方法是指将道路附着系数设为经验常值0.8，进而利用车辆的动力学模型建立状态方程和观测方程，再利用标准扩展卡尔曼滤波递推估计出纵向和侧向车速。

（1）单附着系数路面仿真

道路附着系数设置为0.55，仿真时间为50s，为检验算法在较高机动环境下的估计效果，设置汽车的方向盘转角按正弦规律变化，如附图4所示，纵向车速如附图5所示。表1列出了对于整个过程利用普通扩展卡尔曼滤波算法和本发明方法估计车辆纵向、侧向速度和道路附着系数的统计结果对比，表中的误差均是相对于Carsim输出的相应参考值而言的（如本发明方法的纵向速度误差就表示利用本方明方法估计出的纵向速度相对于Carsim输出的纵向速度参考值的误差）。另外需指出的是，上述两种方法的具体含义如下：普通的扩展卡尔曼滤波方法是指将道路附着系数设为经验常值0.8，进而利用车辆的动力学模型建立状态方程和观测方程，再利用标准扩展卡尔曼滤波递推估计出纵向和侧向速度；本发明方法是指利用本发明提出的基于交互多模型的估计方法来估计车辆纵向和侧向速度的方法。

表1两种方法在单附着系数路面上估计效果的对比表

表中“--”表示普通的扩展卡尔曼滤波方法无法推算的项

图6给出了利用本发明方法所估计出的道路附着系数值，图中估计结果用灰虚线代表，Carsim输出真值用黑实线代表为进一步说明两种估计方法的优劣，以纵向速度为例，图7给出了本发明方法的纵向速度估计误差，图8给出基于扩展卡尔曼滤波算法的纵向速度估计误差。

由表1的对比（尤其是标准差）以及图7～图8，可以看出本发明方法相对于普通扩展卡尔曼滤波方法在纵向车速和侧向速度的估计方面精度有了大幅的提高。另外，根据表1及图6，还可以看出本发明方法能够实时的估计出道路附着系数，并且具有较高的精度，这也是普通扩展卡尔曼滤波方法所无法实现的。

（2）附着系数突变路面仿真

为检验本发明所提出算法对附着系数突变路面的适应性，道路附着系数设置为由0.8到0.4再到0.6跃变，仿真时间为50s，方向盘转角变化如附图9所示，同时所设置的车辆纵向速度在不断地做加速、制动减速和匀速等变化，以检验本发明所提出方法既可适应一般机动环境也可适应较高机动环境下车辆运行状态的准确估计，纵向速度如图10所示。表2列出了对于整个过程利用普通扩展卡尔曼滤波算法和本发明方法估计车辆纵向、侧向速度和道路附着系数的统计结果对比，表中的误差均是相对于Carsim输出的相应参考值而言的（如本发明方法的纵向速度误差就表示利用本方明方法估计出的纵向速度相对于Carsim输出的纵向速度参考值的误差）。

表2两种方法在附着系数突变路面上估计效果的对比表

表中“--”表示普通的扩展卡尔曼滤波方法无法推算的项

图11给出了利用本发明方法所估计出的道路附着系数值,图中估计结果用灰虚线代表，Carsim输出真值用黑实线代表。

由表2的对比（尤其是标准差）以及图11，可以看出本发明方法相对于普通扩展卡尔曼滤波方法在纵向车速和侧向速度的估计方面精度有了大幅的提高。另外，根据表2及图11，还可以看出本发明方法在道路附着系数突变情况下，能够迅速的识别出道路附着系数的改变，并且具有较高的估计精度，这也是普通扩展卡尔曼滤波方法所无法实现的。

综上，即使在道路附着系数突变的环境下，本发明提出的方法能够准确地估计出车辆纵向前进速度、侧向速度，并实时估计出道路附着系数信息，且既可适应一般机动环境也可适应较高机动环境，这些信息可满足有关汽车主动安全控制的需要。

Claims (1)

1.一种车速与道路附着系数的联合估计方法，其特征在于：本方法是针对前轮转向四轮汽车，基于非线性整车动力学模型和轮胎纵向力模型，在不同道路附着系数条件下，分别建立不同的多个卡尔曼滤波模型，同时利用车载轮速和方向盘转角传感器信息来确定建立各卡尔曼滤波系统的外部输入量和观测量。进一步通过交互多模型算法实现不同道路附着系数条件下对车辆纵向、侧向车速的自适应估计，并根据交互多模型算法中计算出的各卡尔曼滤波模型的模型概率实现道路附着系数的实时估计，达到全面自适应的效果；

具体步骤包括：

1）建立扩展卡尔曼滤波的状态方程和观测方程

针对道路附着系数分别为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0，即μ_j＝10×j（j=1,2,...10）时，分别建立10个不同的扩展卡尔曼滤波模型，其中μj为针对于不同模型的道路附着系数；这10个模型具有相同的形式，其区别仅在于道路附着系数具体取值的不同；根据三自由度的汽车非线性动力学模型建立扩展卡尔曼滤波的系统状态方程，第j（j=1,2,…10）个模型离散化后的卡尔曼滤波的状态方程的矩阵形式表示为：

X_j(k)＝f_j(X_j(k-1),U_j(k-1),W_j(k-1),γ_j(k-1))            （1）

式（1）中，下标j表示第j个模型（j=1,2,...10），k表示离散化时刻；这10个模型具有相同的系统状态向量，该系统状态向量为X_j＝[x₁ x₂ x₃]′，其中，x₁＝v_x，x₂＝v_y，x₃＝r，v_x、v_y及r分别是汽车的纵向前进速度、侧向速度和横摆角速度，本发明中上角标'表示对矩阵转置；系统外输入向量为U_j＝[u₁ u₂ u₃]′，其中，u₁＝δ，u₂＝F_{j_xf}，u₃＝F_{j_xr}，δ是前轮转向角，F_{j_xf}表示第j个模型（j=1,2,...10）中作用在单个前轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个前轮上的纵向力，F_{j_xr}是第j个模型（j=1,2,...10）中作用在单个后轮上的纵向力，即当道路附着系数为μ_j（j=1,2,...10）时作用在单个后轮上的纵向力；W_j表示零均值的系统高斯白噪声向量且W_j＝[w₁ w₂ w₃]′，其中w₁、w₂及w₃分别表示三个系统高斯白噪声分量；γ_j表示系统外输入对应的零均值高斯白噪声向量且 ${\mathrm{\γ}}_j={\left[\begin{array}{ccc}{w}_{\mathrm{\δ}}& w_{F_{\mathrm{j\; xf}}}& w_{F_{\mathrm{j\; xr}}}\end{array}\right]}^{\′},$ 其中w_δ表示系统外输入δ对应的零均值高斯白噪声，

和

分别表示外输入F_{j_xf}和F_{j_xr}对应的零均值高斯白噪声，这些白噪声隐含在系统状态函数的系统外输入里面；

非线性的系统状态函数向量为

$f_j(X_j,U_j,W_j,{\mathrm{\γ}}_j)=\left[\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_2}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\\ f_{j\_3}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))\end{array}\right],$

其中，

$\begin{array}{c}{f}_{j\_1}(X_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=v_x(k-1)+\frac{T}{M}[{\mathrm{Mv}}_y(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}\frac{v_y(k-1)+\mathrm{ar}(k-1)}{v_x(k-1)}\mathrm{\δ}(k-1)]+\frac{2T}{M}[F_{j\_\mathrm{xf}}(k-1)+F_{j\_\mathrm{xr}}(k-1)]+w_1\\ f_{j\_2}(X_j(k-1),U_j(k-1),w_j(k-1),{\mathrm{\γ}}_j(k-1))=v_y(k-1)+\frac{T}{M}\{-{\mathrm{Mv}}_x(k-1)r(k-1)+\\ 2C_{\mathrm{\αf}}[\mathrm{\δ}(k-1)+\frac{-v_y(k-1)-\mathrm{ar}(k-1)}{v_x(k-1)}]+2C_{\mathrm{\αr}}\frac{\mathrm{br}(k-1)-v_y(k-1)}{v_x(k-1)}\}+\frac{2T}{M}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_2\\ f_{j\_3}(x_j(k-1),U_j(k-1),W_j(k-1),{\mathrm{\γ}}_j(k-1))=r(k-1)+\frac{T}{I_z}\left\{2a{C}_{\mathrm{\αf}}\right[\mathrm{\δ}(k-1)-\frac{(v_y(k-1)+\mathrm{ar}(k-1))}{v_x(k-1)}]\\ -2b{C}_{\mathrm{ar}}\frac{[\mathrm{br}(k-1)-v_y(k-1)]}{v_x(k-1)}\}+\frac{2\mathrm{aT}}{I_z}{F}_{j\_\mathrm{xf}}(k-1)\mathrm{\δ}(k-1)+w_3\end{array}$

在上述表达式中，M和I_z分别是车辆的质量和绕过质心垂向轴的转动惯量，a是汽车前轮轮轴中心到质心的距离，b是汽车后轮轮轴中心到质心的距离，C_αf、C_αr分别表示前、后轮胎的侧偏刚度，T表示离散的周期，其典型值为10毫秒、20毫秒、50毫秒或100毫秒；W_j对应的系统噪声协方差阵Q_j为：

$Q_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{w_1}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{w_2}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{w_3}}^2\end{array}\right],$ 其中

及分别表示系统高斯白噪声w₁、w₂及w₃对应的方差；γ_j对应的系统外部输入噪声的协方差阵为 ${\mathrm{\Γ}}_j=\left[\begin{array}{ccc}{{\mathrm{\σ}}_{\mathrm{\δ}}}^2& 0& 0\\ 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xf}}}}^2& 0\\ 0& 0& {{\mathrm{\σ}}_{F_{j\_\mathrm{xr}}}}^2\end{array}\right],$

和分别表示w_δ、

和

对应的方差；轮胎纵向力F_{j_xf}和F_{j_xr}根据非线性刷子轮胎模型来确定，轮胎模型中道路附着系数μ_j（j=1,2,...10）取值的不同是10个模型的区别所在；

用s_xq（q=f，r）表示车辆纵向滑移率，即又可分为前轮轴纵向滑移率s_xf和后轮轴纵向滑移率s_xr，下角标q取f或r，f或r分别表示前或后轮轴，s_xq计算方法为：

s_xq＝(ω_qR-v_xq)/max(ω_qR,v_xq)(q＝f,r)               （2）

式（2）中，R表示车轮轮胎半径；v_xf和v_xr分别表示前、后轮轴上沿轮胎方向的速度，v_xf和v_xr可统一记为v_xq(q=f,r)；max表示求最大值；ω_f表示前轮轴上两个车轮的旋转角速度等效折算到前轮轴上的旋转角速度；ω_r表示后轮轴上两个车轮旋转角速度等效折算到后轮轴上的旋转角速度，ω_f和ω_r可统一记为ω_q(q=f,r)且

$\begin{array}{c}{\mathrm{\ω}}_f=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{fR}}+{\mathrm{\ω}}_{\mathrm{fL}})\\ {\mathrm{\ω}}_r=\frac{1}{2}({\mathrm{\ω}}_{\mathrm{rR}}+{\mathrm{\ω}}_{\mathrm{rL}})\end{array}---\left(3\right)$

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

v_xq(q=f,r)可按式（4）确定：

$\begin{array}{c}{v}_{\mathrm{xf}}=v_x\cos \mathrm{\δ}+(v_y+\mathrm{ar})\sin \mathrm{\δ}\\ v_{\mathrm{xr}}=v_x\end{array}---\left(4\right)$

进而，轮胎纵向力F_{j_xf}和F_{j_xr}可通过式（5）来确定

$F_{j\_\mathrm{xq}}=C_{\mathrm{xq}}{s}_{\mathrm{xq}}-\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^2}{3{\mathrm{\μ}}_j{F}_{\mathrm{zq}}}+\frac{{\left(C_{\mathrm{xq}}{s}_{\mathrm{xq}}\right)}^3}{27{\left({\mathrm{\μ}}_j{F}_{\mathrm{zq}}\right)}^2}(q=f,\mathrm{rj}=\mathrm{1,2},...10)---\left(5\right)$

式（5）中，C_xf和C_xr分别表示单个前、后轮胎的纵向刚度，统一记为C_xq(q=f,r)；μ_j（j＝1,2...10）表示轮胎和地面间的道路摩擦系数，所建10个模型的区别仅在于其取值的不同，其中，μ₁＝0.1，μ₂＝0.2...μ₁₀＝1.0；

F_zq(q=f,r)表示分配到单个前或后轮上的垂向载荷且可按下式计算

$F_{\mathrm{zf}}=\frac{\mathrm{Mgb}}{(a+b)},F_{\mathrm{zr}}=\frac{\mathrm{Mga}}{(a+b)}---\left(6\right)$

式（6）中，g表示重力加速度；

车辆纵向前进速度和横摆角速度与两个非转向后轮的速度存在以下关系

$\begin{array}{c}{v}_x=(V_{\mathrm{RL}}+V_{\mathrm{RR}})/2\\ r=(V_{\mathrm{RL}}-V_{\mathrm{RR}})/T_W\end{array}---\left(7\right)$

式（7）中，T_W表示后轮轴上两个后轮间的轮距，V_RL和V_RR分别表示左后轮和右后轮的线速度；

第j（j=1,2,...10）个模型的卡尔曼滤波的观测方程的离散化矩阵形式为：

Z_j(k)=H_j(k)X_j(k)+V_j(k)               （8）

式（8）中，Z_j为观测向量，H_j为观测阵，V_j表示与W_j互不相关的零均值观测白噪声向量，且 $Z_j\left(k\right)=\left[\begin{array}{c}{v}_{x\_m}\left(k\right)\\ {\mathrm{\ω}}_{z\_m}\left(k\right)\end{array}\right],H_j\left(k\right)=\begin{array}{}\end{array}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right],V_j=\left[\begin{array}{c}{n}_{v_x}\\ n_{{\mathrm{\ω}}_z}\end{array}\right],$ 其中v_{x_m}(k)和ω_{z_m}(k)分别为通过轮速传感器测量获得的车辆纵向前进速度和横摆角速度；

表示通过轮速传感器测量获得的车辆纵向前进速度的观测噪声且

是均值为0、方差为的高斯白噪声，

表示通过轮速传感器测量获得的横摆角速度的观测噪声且

是均值为0、方差为

的高斯白噪声；V_j对应的观测噪声方差阵R_j可表示为 $R_j=\left[\begin{array}{cc}{{\mathrm{\σ}}_{v_x}}^2& 0\\ 0& {{\mathrm{\σ}}_{{\mathrm{\ω}}_z}}^2\end{array}\right];$

对于式（8）中的测量值v_{x_m}(k)和ω_{z_m}(k)，它们是利用后轮轴上两个轮速传感器测得的角速度乘以轮胎半径得到V_{RL_m}＝R·ω_rL和V_{RR_m}＝R·ω_rR，V_{RL_m}和V_{RR_m}分别表示V_RL和V_RR含有噪声的测量值，进而利用式（7）获得的，即v_{x_m}和ω_{z_m}分别表示v_x和r的含有噪声的测量值且

2）交互多模型估计方法

对于式（1）描述的系统状态方程和式（8）描述的测量方程，可运用交互多模型滤波理论，建立起滤波递推估计过程。具体估计步骤如下：

①交互估计计算

上述十个扩展卡尔曼系统模型之间的转移概率为p_ij，下标i、j（i=1,2...10,j=1,2,3...10）表示从状态i转移到状态j的概率；

则预测第j(j=1,2...10)个模型的模型概率ρ_j(k,k-1):

${\mathrm{\ρ}}_j=(k,k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}{\mathrm{\ρ}}_i(k-1)$

预测混合概率ρ_i|j(k-1):

ρ_i|j(k-1)＝p_ijρ_i(k-1)/ρ_j(k,k-1)

则交互估计后第j个滤波器的输入为：

$X_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{X}_i(k-1){\mathrm{\ρ}}_{i|j}(k-1)$

$P_{0j}(k-1)=\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{i|j}(k-1)\{P_i(k-1)+[X_i(k-1)-X_{0j}(k-1)\left]\right[X_i(k-1)-X_{0j}(k-1){]}^{\′}\}$

②每个模型滤波器对于式（1）和式（8）所描述的状态方程和观测方程，运用扩展卡尔曼滤波理论，各自进行标准扩展卡尔曼滤波递推，该递推过程包括时间更新和测量更新，第j(j=1,2,...10)个模型的滤波过程如下：

时间更新：

状态一步预测方程X_j(k,k-1)＝f_j(X_0j(k-1),U_j(k-1),0,0)

一步预测误差方差阵：

P_j(k,k-1)＝A_j(k-1)P_0j(k-1)(A_j(k-1))′+B_j(k-1)Γ_j(k-1)(B_j(k-1))′+Q_j(k-1)

其中，A_j、B_j分别是系统状态函数向量f_j对状态向量X_j和外部输入向量U_j求偏导数的雅可比矩阵，即矩阵A_j和B_j的第m行第n列元素A_{j_[m,n]}和B_{j_[m,n]}可分别通过下式求得：

$A_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂x_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

$B_{j\_[m,n]}=\frac{\∂f_{j\_m}}{\∂u_n}(X_j(k,k-1),U_j(k-1),\mathrm{0,0})(m=\mathrm{1,2,3}n=\mathrm{1,2,3})$

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

$A_{j\_\left[\mathrm{1,1}\right]}=1+T\left[\frac{-2C_{\mathrm{\αf}}(v_y+\mathrm{ar})}{M{v_x}^2}\mathrm{\δ}\right]A_{j\_\left[\mathrm{1,2}\right]}=T[r+\frac{{2C}_{\mathrm{\αf}}}{M_{v_x}}\mathrm{\δ}]$ $A_{j\_\left[\mathrm{1,3}\right]}=T(v_y+\frac{2C_{\mathrm{\αf}}a}{M_{v_x}}\mathrm{\δ})A_{j\_\left[\mathrm{2,1}\right]}=T[-r-2C_{\mathrm{\αr}}\frac{\mathrm{br}-v_y}{M{v_x}^2}+2C_{\mathrm{\αf}}\frac{v_y+\mathrm{ar}}{M{v_x}^2}]$

$A_{j\_\left[\mathrm{2,2}\right]}=1-\frac{2T(C_{\mathrm{\αr}}+C_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}{A}_{j\_\left[\mathrm{2,3}\right]}=T[-v_x+2\frac{(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{{\mathrm{Mv}}_x}]$

$A_{j\_\left[\mathrm{3,1}\right]}=\frac{2T[a{C}_{\mathrm{\αf}}(v_y+\mathrm{ar})+b{C}_{\mathrm{\αr}}(\mathrm{br}-v_y)]}{I_z{v_x}^2}$

$A_{j\_\left[\mathrm{3,2}\right]}=\frac{2T(b{C}_{\mathrm{\αr}}-a{C}_{\mathrm{\αf}})}{I_z{v}_x}{A}_{j\_\left[\mathrm{3,3}\right]}=1-\frac{2T(a^2{C}_{\mathrm{\αf}}+b^2{C}_{\mathrm{\αr}})}{I_z{v}_x}$

$B_{j\_\left[\mathrm{1,1}\right]}=2T{C}_{\mathrm{\αf}}\frac{(v_y+\mathrm{ar})}{M{v}_x}{B}_{j\_\left[\mathrm{1,2}\right]}=\frac{2T}{M}{B}_{j\_\left[\mathrm{1,3}\right]}=\frac{2T}{M}$

$B_{j\_\left[\mathrm{2,1}\right]}=\frac{2T{F}_{j\_\mathrm{xf}}}{M}+\frac{2T{C}_{\mathrm{\αf}}}{M}{B}_{j\_\left[\mathrm{2,2}\right]}=\frac{2\mathrm{T\δ}}{M}{B}_{j\_\left[\mathrm{2,3}\right]}=0$

$B_{j\_\left[\mathrm{3,1}\right]}=\frac{2\mathrm{Ta}}{I_z}{C}_{\mathrm{\αf}}+\frac{2\mathrm{Ta}}{I_z}{F}_{j\_\mathrm{xf}}{B}_{j\_\left[\mathrm{3,2}\right]}=\frac{2\mathrm{Ta\δ}}{I_z}{B}_{j\_\left[\mathrm{3,3}\right]}=0$

测量更新：

滤波增益矩阵：K_j(k)＝P_j(k,k-1)(H_j(k))′(S_j(k))^-1

S_j(k)＝H_j(k)P_j(k,k-1)(H_j(k))′+R_j(k)

状态估计：X_j(k)＝X_j(k,k-1)+K_j(k)(Z_j(k)-H_j(k)X_j(k,k-1))

估计误差方差阵：P_j(k)＝P_j(k,k-1)-K_j(k)S_j(k)(K_j(k))′

③模型概率更新

在每个模型完成上一步的更新之后，利用最大似然函数Λ_j(k)计算新的模型概率ρ_j(k)，最大似然函数计算如下：

${\mathrm{\Λ}}_j\left(k\right)=\frac{\exp \{-\frac{1}{2}{(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))}^{\′}{\left(S_j\left(k\right)\right)}^{-1}(Z_j\left(k\right)-H_j\left(k\right)X_j(k,k-1))\}}{|2\mathrm{\π}{S}_j\left(k\right){|}^{-\frac{1}{2}}}$

因此，模型j在k时刻的模型概率由贝叶斯定理给出：

${\mathrm{\ρ}}_j\left(k\right)=\frac{{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_j(k,k-1)}{\underset{i=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\Λ}}_j\left(k\right){\mathrm{\ρ}}_i(k,k-1)}$

④估计组合

在计算出各模型为正确的后验概率之后，对所有滤波器的状态估计进行概率加权并求和，权系数为模型正确的后验概率，得到最终的状态估计为：

$X\left(k\right)=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{X}_j\left(k\right){\mathrm{\ρ}}_j\left(k\right),$ 其中， $X\left(k\right)={\left[\begin{array}{ccc}\stackrel{\‾}{v_x}& \stackrel{\‾}{v_y}& \stackrel{\‾}{r}\end{array}\right]}^{\′},$ 各状态变量的上标“-”表示各状态量的最终滤波估计值，即X(k)内各状态变量依次分别表示估计组合后的纵向车速、侧向车速和横摆角速度；

同时，由于各模型的区别在于各模型所设定的道路附着系数的具体取值不同，即各模型的μ_j的取值不同，因此，对各模型所设定的附着系数进行概率加权即可得出最终滤波估计出的当前的道路附着系数μ：

$\mathrm{\μ}=\underset{j=1}{\overset{10}{\mathrm{\Σ}}}{\mathrm{\μ}}_j{\mathrm{\ρ}}_j\left(k\right)---\left(9\right).$

Images