本发明属于汽车技术领域，具体的说是一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法。包括以下步骤：步骤一、基于车辆动力学模型对车辆的状态参数进行求解；步骤二、基于卡尔曼滤波估计出车辆行驶过程中轮胎受到的轮胎力；步骤三、采用最小二乘法对路面附着系数进行了估计；步骤四、通过斜率法对步骤三中低滑移率下的路面附着系数估计结果进行修正。本发明有效地降低了轮胎力观测值的不确定性，能够精确的估计在车辆滑移率较高时的路面附着系数，在车辆滑移率较低时对估计的路面附着系数结果进行合理的修正。
1.一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，其特征在于，包括以下步骤：步骤一、基于车辆动力学模型对车辆的状态参数进行求解；
步骤二、基于卡尔曼滤波估计出车辆行驶过程中轮胎受到的轮胎力；
步骤四、通过斜率法对步骤三中低滑移率下的路面附着系数估计结果进行修正；
步骤一中所述状态参数包括车速、制动压力、前轮转角、轮胎滑移率、侧偏角和车轮垂向力；所述车速、制动压力、前轮转角由传感器获得；所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算出；
所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算的具体方法如下：
11)建立包含车辆纵向、侧向和横摆运动的三自由度模型，得车辆动力学模型方程为：max＝[cosδ cosδ 1 1 ‑sinδ 0]Ftire (1)may＝[sinδ sinδ 0 0 cosδ 1]Ftire (2)T
式中，Ftire＝[Fxfl Fxfr Fxrl Fxrr FyfFyr]表示轮胎力集合；Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力；Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyrl表示后左轮侧向力；Fyrr表示后右轮侧向力；m表示整车质量，ax表示车辆纵向加速度；ay表示车辆侧向加速度；δ表示车辆前轮转角；Iz表示车辆绕z轴的转动惯量；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；tw表示轮距；表示车辆质心处横摆角加速度；
式中，i＝fl，fr，rl，rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Iwi表示轮胎旋转惯性质量；表示车轮转动角加速度；Tdi表示车轮的驱动转矩；Tbi表示车轮的驱动转矩与制动转矩；Fxi表示车轮受到地面的纵向反作用力；r表示车辆滚动半径；TFi＝(a+bvi)Fzir表示车轮受到的滚动阻力矩；a和b表示滚动阻力系数参数值；vi表示车轮中心处的纵向速度；Fzi表示车轮受到地面的垂向反作用力；
式中，Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Te表示整车动力电机输出转矩；i0表示传动系的传动比；ηt表示传动系的传动效率；
14)车辆制动过程中由制动力形成的制动力矩表示为：
式中，i＝fl，fr，rl，rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Tbi表示车轮的制动转矩；kbi表示制动效能因数；Pwi表示轮缸制动压力；
15)车轮受到地面的垂向反作用力与车辆加速和减速形成的惯性力存在关系，表示为：式中，Fzfl表示前左车轮受到地面的垂向反作用力；Fzfr表示前右车轮受到地面的垂向反作用力；Fzrl表示后左车轮受到地面的垂向反作用力；Fzrr表示后右车轮受到地面的垂向反作用力；m表示整车质量；ax表示车辆纵向加速度；ay表示车辆侧向加速度；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；hg表示质心高度；tw表示轮距；
16)选用前左轮纵向力，前右轮纵向力，后左轮纵向力，后右轮纵向力，前轮总侧向力，后轮总侧向力，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统状态变量，即选用车辆纵向加
速度，车辆侧向加速度，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统输出，那么得到车辆动力学空间方程为：
式中，X表示系统状态变量；表示系统状态变量的微分；u表示系统输入；Φ表示系统状态变量矩阵；B表示系统输入矩阵；Y表示系统输出；H表示系统输出矩阵，且有：Hx＝[cosδ cosδ 1 1 ‑sinδ 0]，
Hy＝[sinδ sinδ 0 0 cosδ 1]，
式中，Φ表示系统状态变量矩阵；B表示系统输入矩阵；H表示系统输出矩阵；r表示车辆滚动半径；I表示单位矩阵；Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；
Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Tbfl表示前左车轮的制动转矩；Tbfr表示前右车轮的制动转矩；Tbrl表示后左车轮的制动转矩；Tbrr表示后右车轮的制动转矩；Iwfl表示前左车轮的轮胎旋转惯性质量；Iwfr表示前右车轮的轮胎旋转惯性质量；Iwrl表示后左车轮的轮胎旋转惯性质量；Iwrr表示后右车轮的轮胎旋转惯性质量；Iz表示车辆绕Z轴的转动惯量；Hx表示纵向力矩阵；Hy表示侧向力矩阵； 表示横摆角矩阵；lf表示车辆质心至前轴的距离；lr表示车辆质心至后轴的距离；δ表示车辆前轮转角；tw表示轮距；m表示整车质量；TFfl表示前左车轮受到的滚动阻力矩；TFfr表示前右车轮受到的滚动阻力矩；TFrl表示后左车轮受到的滚动阻力矩；TFrr表示后右车轮受到的滚动阻力矩；
17)采用刷子模型表征车辆轮胎的纵向力和车辆轮胎的侧向力之间的联系，表示为：式中，i＝fl，fr，rl，rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Fxi表示车轮受到地面的纵向反作用力；Fyi表示车轮受到地面的侧向反作用力；Cxi表示轮胎的纵向刚度；Cαi表示轮胎的侧向刚度；κi表示轮胎的滑移率；αi表示轮胎侧偏角；fi表示刷子模型刚度系数；Fti表示刷子模型表征力；μ表示路面附着系数；Fzi表示车轮受到地面的垂向反作用力；
18)根据车轮旋转角速度、前轮转角、车速以及车辆结构参数，计算各个轮胎的滑移率：式中，κfl表示前左轮滑移率；κfr表示前右轮滑移率；κrl表示后左轮滑移率；κrr表示后右轮滑移率；ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度；r表示车辆滚动半径；vx表示车辆纵向速度；vy表示车辆侧向速度；tw表示轮距；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；lf表示车辆质心到前轴距离；表示车辆质心处横摆角速度；
19)智能电动汽车在行驶过程中，对应轮胎侧偏角为：
式中，αfl表示前左轮轮胎侧偏角；αfr表示前右轮轮胎侧偏角；αrl表示后左轮轮胎侧偏角；αrr表示后右轮轮胎侧偏角；表示车辆质心处横摆角速度；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；β表示车辆质心侧偏角；vx表示车辆纵向速度；
2.根据权利要求1所述的一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，其特征在于，所述步骤二的具体方法如下：
21)将连续的车辆动力学空间方程(8)转换为离散系统状态空间方程，即：式中，X(k)表示第k时刻离散系统的状态变量；u(k)表示第k时刻离散系统的输入；X(k+
1)表示第k+1时刻离散系统的状态变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；Y(k)表示第k时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵；w(k)表示第k时刻离散系统的激励噪声；v(k)表示第k时刻离散系统的观测噪声；
22)假设离散系统的激励噪声w(k)和观测噪声v(k)是均值为零、方差各为Q和R的不相关白噪声，离散系统的初始状态X(0)不相关于激励噪声w(k)和观测噪声v(k)；那么卡尔曼滤波的一步预测为：式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值； 表示第k时刻的离散系统变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；u(k)表示第k时刻离散系统的输入；
式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值；K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；Y(k+1)表示第k+1时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵； 表示经校正后第k+1时刻卡尔曼滤波的估计值；
式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Hk表示第k时刻离散系统的输出矩阵；R表示观测噪声v(k)的方差；
式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第k时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；
26)在一步预测协方差矩阵中引入遗忘因子进行修正，即：
式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第k时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；s表示遗忘因子；
P(k+1|k+1)＝(I‑K(k+1)Hk)P(k+1|k) (24)式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；Hk表示第k时刻离散系统的输出矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；P(k+1|k+1)为经校正后的协方差矩阵；I表示单位矩阵；
2 7) 通 过 式 ( 1 9) ‑ ( 2 4) 迭 代 ，估 算 获 得 离 散 系 统的 状 态 变 量中的轮胎纵向力和侧向力；
其中，Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力；Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyrl表示后左轮侧向力；Fyrr表示后右轮侧向力；
ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度，表示车辆质心处横摆角速度；
3.根据权利要求1所述的一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，其特征在于，所述步骤三的具体方法如下：
式中，y(k)＝[Fx，Fy]表示第k时刻的轮胎力；Fx表示车轮受到地面的纵向反作用力；Fy表示车轮受到地面的侧向反作用力；f(k，θ(k))表示第k时刻的刷子轮胎模型；θ(k)＝[Cx，Cα，Tμ] 表示第k时刻刷子轮胎模型的状态向量；Cx表示轮胎的纵向刚度；Cα表示轮胎的侧向刚度；μ表示路面附着系数；
32)对y(k)进行一阶泰勒展开，并忽略高阶项，得：
y(k)≈F(k)(θ(k)‑θ(k‑1))+f(θ(k‑1)，k‑1) (26)式中，y(k)表示第k时刻的轮胎力；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；θ(k)为第k时刻刷子轮胎模型的状态向量；f(θ(k‑1)，k‑1)表示第k‑1时刻的刷子轮胎模型；F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置，且有：式中，F(k)表示y(k)对θ(k)的雅可比矩阵的转置；θ(k)为第k时刻刷子轮胎模型的状态向量；f(k，θ)表示刷子轮胎模型；θ表示刷子轮胎模型的状态向量；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；
令z(k)＝y(k)‑f(θ(k‑1)，k‑1)+F(k)θ(k‑1)，那么：z(k)≈F(k)θ(k) (28)
式中，F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k时刻的轮胎力；
θ(k)第k时刻刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；
33)通过最小二乘法对路面附着系数进行估计，设定最小二乘法代价函数为：式中， 表示基于刷子轮胎模型的递归最小二乘法代价函数；Λ表示最小二乘法的遗忘因子矩阵；z(i)表示第i时刻可测量轮胎力矩阵；F(i)表示第i时刻y(i)对θ(i)的雅可比矩阵的转置；θ(i)为第i时刻刷子轮胎模型的状态向量； 表示第i时刻估计的刷子轮胎模型的状态向量；
34)为了使得代价函数最小，基于刷子轮胎模型的递归最小二乘法为：
L(k)＝P(k‑1)F(k)(I+F(k)P(k‑1)F(k)) (31)‑1 ‑1
P(k)＝Λ (I‑L(k)F(k))P(k‑1)Λ (32)式中， 表示第k时刻估计的刷子轮胎模型的状态向量； 表示第k‑1时刻估计的刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；F(k‑1)表示第k‑1时刻y(k‑1)对θ(k‑1)的雅可比矩阵的转置；y(k‑1)表示第k‑1时刻的轮胎力，θ(k‑1)第k‑1时刻刷子轮胎模型的状态向量，L(k)表示第k时刻基于刷子模型的递推最小二成法增益矩阵；F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k时刻的轮胎力；P(k)表示第k时刻基于刷子模型的最小二乘法更新矩阵；P(k‑1)表示第k‑1时刻基于刷子模型的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵，I表示单位矩阵；
35)通过(30)‑(32)迭代，估算获得估计的刷子轮胎模型的状态向量中的路面附着系数；其中， 表示基于刷子模型估计的轮胎的纵向刚度， 表示基于刷子模型估计的轮胎的侧向刚度，表示基于刷子模型估计的路面附着系数。
4.根据权利要求1所述的一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，其特征在于，所述步骤四的具体方法如下：
41)在滑移率小时，附着系数μ与滑移率κ近似为正比例关系，即：
式中，κ (k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；μs(k)表示第k时刻用于斜率法真实的路面附着系数；
式中， 表示基于斜率法的递归最小二乘法代价函数；Λ表示最小二乘法的遗忘‑1因子矩阵；κ (i)表示第i时刻轮胎的滑移率的倒数；Fz(i)表示第i时刻车轮受到地面的垂向反作用力，Fx(i)表示第i时刻车轮受到地面的纵向反作用力， 为第i时刻斜率法估算的路面附着系数；
43)为了使得代价函数最小，设计基于斜率法的递归最小二乘法为：
式中， 表示第k时刻斜率法估算的路面附着系数； 表示第k‑1时刻斜率法‑1
估算的路面附着系数；κ (k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；Ls(k)表示第k时刻基于斜率法的递推最小二成法增益矩阵；Ps(k)表示第k时刻基于斜率法的最小二乘法更新矩阵；Ps(k‑1)表示第k‑1时刻基于斜率法的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵；I表示单位矩阵；
44)通过式(35)‑(37)迭代，估算获得斜率法估算的路面附着系数最终，修正后的路面附着系数估计结果表示为：
式中， 表示修正后的路面附着系数估计结果；κ表示轮胎滑移率； 表示基于斜率法估计的路面附着系数；表示基于刷子模型估计的路面附着系数。
一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计
方法
技术领域
[0001] 本发明属于汽车技术领域，具体的说是一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法。
背景技术
[0002] 智能电动汽车在地面上行驶一方面受到车辆底层驱动系统的驱动力和制动系统的制动力限制，另一方面也受到不同路面附着条件的约束。路面附着系数是车—路系统中最为关键的参数之一，它对智能电动汽车的行驶安全性能有着至关重要的影响。
[0003] 依托车载摄像头、激光雷达、超声波雷达等是近年来进行路面识别的一种新的方法，但是这些设备成本较高，实际应用不确定性较强，因此并未得到广泛的推广。依托已有车载传感器信息，采用观测估计理论对路面附着系数进行准确辨识是一种更接近实际低成本应用的手段。然而实际应用中，应用于路面附着系数估计的轮胎力不确定性较强，需要采用合适的策略对其进行合理的预测估计，避免路面附着系数估计发散。此外，在轮胎滑移率过低时，基于经典轮胎模型估计的路面附着系数可靠性较低，对车辆运动控制造成安全隐患。
发明内容
[0004] 本发明提供了一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，该方法有效地降低了轮胎力观测值的不确定性，能够精确的估计在车辆滑移率较高时的路面附着系数，在车辆滑移率较低时对估计的路面附着系数结果进行合理的修正。
[0005] 本发明技术方案结合附图说明如下：
[0006] 一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，包括以下步骤：
[0007] 步骤一、基于车辆动力学模型对车辆的状态参数进行求解；
[0008] 步骤二、基于卡尔曼滤波估计出车辆行驶过程中轮胎受到的轮胎力；
[0009] 步骤三、采用最小二乘法对路面附着系数进行了估计；
[0010] 步骤四、通过斜率法对步骤三中低滑移率下的路面附着系数估计结果进行修正。
[0011] 步骤一中所述状态参数包括车速、制动压力、前轮转角、轮胎滑移率、侧偏角和车轮垂向力；所述车速、制动压力、前轮转角由传感器获得；所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算出。
[0012] 所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算的具体方法如下：
[0013] 11)建立包含车辆纵向、侧向和横摆运动的三自由度模型，得车辆动力学模型方程为：
[0014] max＝[cosδcosδ1 1‑sinδ0]Ftire (1)
[0015] may＝[sinδsinδ0 0 cosδ1]Ftire (2)
[0016]
[0017] 式中，Ftire＝[Fxfl Fxfr Fxrl Fxrr Fyf Fyr]T表示轮胎力集合；Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力；Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyrl表示后左轮侧向力；Fyrr表示后右轮侧向力；m表示整车质量，ax表示车辆纵向加速度；ay表示车辆侧向加速度；δ表示车辆前轮转角；Iz表示车辆绕z轴的转动惯量；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；tw表示轮距； 表示车辆质心处横摆角加速度；
[0018] 12)根据动力学平衡关系得到车轮的动力学方程为：
[0019]
[0020] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Iwi表示轮胎旋转惯性质量； 表示车轮转动角加速度；Tdi表示车轮的驱动转矩；Tbi表示车轮的驱动转矩与制动转矩；Fxi表示车轮受到地面的纵向反作用力；r表示车辆滚动半径；TFi＝(a+bvi)Fzir表示车轮受到的滚动阻力矩；a和b表示滚动阻力系数参数值；vi表示车轮中心处的纵向速度；Fzi表示车轮受到地面的垂向反作用力；
[0021] 13)智能电动汽车的四轮驱动力矩可以表示为：
[0022]
[0023] 式中，Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Te表示整车动力电机输出转矩；i0表示传动系的传动比；ηt表示传动系的传动效率；
[0024] 14)车辆制动过程中由制动力形成的制动力矩表示为：
[0025] Tbi＝kbiPwi (6)
[0026] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Tbi表示车轮的制动转矩；kbi表示制动效能因数；Pwi表示轮缸制动压力；
[0027] 15)车轮受到地面的垂向反作用力与车辆加速和减速形成的惯性力存在关系，表示为：
[0028]
[0029] 式中，Fzfl表示前左车轮受到地面的垂向反作用力；Fzfr表示前右车轮受到地面的垂向反作用力；Fzrl表示后左车轮受到地面的垂向反作用力；Fzrr表示后右车轮受到地面的垂向反作用力；m表示整车质量；ax表示车辆纵向加速度； ay表示车辆侧向加速度；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；hg表示质心高度；tw表示轮距；
[0030] 16)选用前左轮纵向力，前右轮纵向力，后左轮纵向力，后右轮纵向力，前轮总侧向力，后轮总侧向力，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统状态变量，即选用车辆纵向加
速度，车辆侧向加速度，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统输出，
那么得到车辆动力学空间方程为：
[0031]
[0032] 式中，X表示系统状态变量； 表示系统状态变量的微分；u表示系统输入；Φ表示系统状态变量矩阵；B表示系统输入矩阵；Y表示系统输出；H表示系统输出矩阵，且有：
[0033]
[0034]
[0035] Hx＝[cosδcosδ1 1 ‑sinδ0],
[0036] Hy＝[sinδsinδ0 0 cosδ1],
[0037]
[0038] 式中，Φ表示系统状态变量矩阵；B表示系统输入矩阵；H表示系统输出矩阵；r表示车辆滚动半径；I表示单位矩阵；Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Tbfl表示前左车轮的制动转矩；Tbfr表示前右车轮的制动转矩；Tbrl表示后左车轮的制动转矩；Tbrr表示后右车轮的制动转矩；Iwfl表示前左车轮的轮胎旋转惯性质量；Iwfr表示前右车轮的轮胎旋转惯性质量；Iwrl表示后左车轮的轮胎旋转惯性质量；Iwrr表示后右车轮的轮胎旋转惯性质量；Iz表示车辆绕Z轴的转动惯量；Hx表示纵向力矩阵；Hy表示侧向力矩阵； 表示横摆角矩阵；lf表示车辆质心至前轴的距离；lr表示车辆质心至后轴的距离；δ表示车辆前轮转角；tw表示轮距；m表示整车质量；TFfl表示前左车轮受到的滚动阻力矩；TFfr表示前右车轮受到的滚动阻力矩；
TFrl表示后左车轮受到的滚动阻力矩；TFrr表示后右车轮受到的滚动阻力矩；
[0039] 17)采用刷子模型表征车辆轮胎的纵向力和车辆轮胎的侧向力之间的联系，表示为：
[0040]
[0041]
[0042]
[0043]
[0044] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Fxi表示车轮受到地面的纵向反作用力；Fyi表示车轮受到地面的侧向反作用力；Cxi表示轮胎的纵向刚度；Cαi表示轮胎的侧向刚度；κi表示轮胎的滑移率；αi表示轮胎侧偏角；fi表示刷子模型刚度系数；Fti表示刷子模型表征力；μ表示路面附着系数；Fzi表示车轮受到地面的垂向反作用力；
[0045] 18)根据车轮旋转角速度、前轮转角、车速以及车辆结构参数，计算各个轮胎的滑移率：
[0046]
[0047]
[0048]
[0049]
[0050] 式中，κfl表示前左轮滑移率；κfr表示前右轮滑移率；κrl表示后左轮滑移率；κrr表示后右轮滑移率；ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度；r表示车辆滚动半径；vx表示车辆纵向速度；vy表示车辆侧向速度；tw表示轮距；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；lf表示车辆质心到前轴距离； 表示车辆质心处横摆角速度；
[0051] 19)智能电动汽车在行驶过程中，对应轮胎侧偏角为：
[0052]
[0053] 式中，αfl表示前左轮轮胎侧偏角；αfr表示前右轮轮胎侧偏角；αrl表示后左轮轮胎侧偏角；αrr表示后右轮轮胎侧偏角； 表示车辆质心处横摆角速度； lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；β表示车辆质心侧偏角；vx表示车辆纵向速度；
[0054] 所述步骤二的具体方法如下：
[0055] 21)将连续的车辆动力学空间方程(8)转换为离散系统状态空间方程，即：
[0056]
[0057] 式中，X(k)表示第k时刻离散系统的状态变量；u(k)表示第k时刻离散系统的输入；X(k+1)表示第k+1时刻离散系统的状态变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；Y(k)表示第k 时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵；w(k)表示第k时刻离散系统的激励噪声；v(k)表示第k时刻离散系统的观测噪声；
[0058] 22)假设离散系统的激励噪声w(k)和观测噪声v(k)是均值为零、方差各为 Q和R的不相关白噪声，离散系统的初始状态X(0)不相关于激励噪声w(k)和观测噪声v(k)；那么卡尔曼滤波的一步预测为：
[0059]
[0060] 式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值； 表示第k时刻的离散系统变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；u(k) 表示第k时刻离散系统的输入；
[0061] 23)离散系统状态更新矩阵为：
[0062]
[0063] 式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值；K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵； Y(k+1)表示第k+1时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵； 表示经校正后第k+1
时刻卡尔曼滤波的估计值；
[0064] 24)卡尔曼滤波增益矩阵表示为：
[0065]
[0066] 式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Hk表示第k 时刻离散系统的输出矩阵；R表示观测噪声v(k)的方差；
[0067] 25)一步预测协方差矩阵为：
[0068]
[0069] 式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第 k时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；
[0070] 26)在一步预测协方差矩阵中引入遗忘因子进行修正，即：
[0071]
[0072] 式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第k 时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；s表示遗忘因子；
[0073] 协方差矩阵更新方程为：
[0074] P(k+1|k+1)＝(I‑K(k+1)Hk)P(k+1|k) (24)
[0075] 式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；Hk表示第k时刻离散系统的输出矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对 k+1时刻的预测协方差矩阵；P(k+1|k+1)为经校正后的协方差矩阵；I表示单位矩阵；
[0076] 27)通过式 (19)‑ (2 4) 迭代 ，估算获得 离散系统的状态变 量中的轮胎纵向力和
侧向力；
[0077] 其中，Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力； Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyrl表示后左轮侧向力；Fyrr表示后右轮侧向力；ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度，表示车辆质心处横摆角速度；
[0078] 所述步骤三的具体方法如下：
[0079] 31)对刷子模型进行改写，即：
[0080] y(k)＝f(k,θ(k)) (25)
[0081] 式中，y(k)＝[Fx,Fy]表示第k时刻的轮胎力；Fx表示车轮受到地面的纵向反作用力；Fy表示车轮受到地面的侧向反作用力；f(k,θ(k))表示第k时刻的刷子轮胎模型；θ(k)＝T[Cx,Cα,μ]表示第k时刻刷子轮胎模型的状态向量；Cx表示轮胎的纵向刚度；Cα表示轮胎的侧向刚度；μ表示路面附着系数；
[0082] 32)对y(k)进行一阶泰勒展开，并忽略高阶项，得：
[0083] y(k)≈F(k)(θ(k)‑θ(k‑1))+f(θ(k‑1),k‑1) (26)
[0084] 式中，y(k)表示第k时刻的轮胎力；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；θ(k)为第k时刻刷子轮胎模型的状态向量；f(θ(k‑1),k‑1)表示第k‑1时刻的刷子轮胎模型；
F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置，且有：
[0085]
[0086] 式中，F(k)表示y(k)对θ(k)的雅可比矩阵的转置；θ(k)为第k时刻刷子轮胎模型的状态向量；f(k,θ)表示刷子轮胎模型；θ表示刷子轮胎模型的状态向量；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；
[0087] 令z(k)＝y(k)‑f(θ(k‑1),k‑1)+F(k)θ(k‑1)，那么：
[0088] z(k)≈F(k)θ(k) (28)
[0089] 式中，F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k 时刻的轮胎力；θ(k)第k时刻刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；
[0090] 33)通过最小二乘法对路面附着系数进行估计，设定最小二乘法代价函数为：
[0091]
[0092] 式中， 表示基于刷子轮胎模型的递归最小二乘法代价函数；Λ表示最小二乘法的遗忘因子矩阵；z(i)表示第i时刻可测量轮胎力矩阵；F(i)表示第i时刻y(i)对θ(i)的雅可比矩阵的转置；θ(i)为第i时刻刷子轮胎模型的状态向量； 表示第i时刻估计的刷子轮胎模型的状态向量；
[0093] 34)为了使得代价函数最小，基于刷子轮胎模型的递归最小二乘法为：
[0094]
[0095] L(k)＝P(k‑1)FT(k)(I+F(k)P(k‑1)FT(k))‑1 (31)
[0096] P(k)＝Λ‑1(I‑L(k)F(k))P(k‑1)Λ‑1 (32)
[0097] 式中， 表示第k时刻估计的刷子轮胎模型的状态向量； 表示第k‑1时刻估计的刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；F(k‑1)表示第k‑1时刻y(k‑1)对θ(k‑1)的雅可比矩阵的转置；y(k‑1) 表示第k‑1时刻的轮胎力，θ(k‑1)第k‑1时刻刷子轮胎模型的状态向量，L(k) 表示第k时刻基于刷子模型的递推最小二成法增益矩阵；F(k)表示第k时刻 y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k时刻的轮胎力；P(k)表示第k 时刻基于刷子模型的最小二乘法更新矩阵；P(k‑1)表示第k‑1时刻基于刷子模型的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵，I表示单位矩阵；
[0098] 35)通过(30)‑(32)迭代，估算获得估计的刷子轮胎模型的状态向量中的路面附着系数；其中， 表示基于刷子模型估计的轮胎的纵向刚度，
表示基于刷子模型估计的轮胎的侧向刚度，表示基于刷子模型估计的路面附着系数。
[0099] 所述步骤四的具体方法如下：
[0100] 41)在滑移率小时，附着系数μ与滑移率κ近似为正比例关系，即：
[0101]
[0102] 式中，κ‑1(k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；μs(k)表示第k时刻用于斜率法真实的路面附着系数；
[0103] 42)设定基于斜率法的最小二乘法代价函数为：
[0104]
[0105] 式中， 表示基于斜率法的递归最小二乘法代价函数；Λ表示最小二乘法的‑1遗忘因子矩阵；κ (i)表示第i时刻轮胎的滑移率的倒数；Fz(i)表示第i时刻车轮受到地面的垂向反作用力，Fx(i)表示第i时刻车轮受到地面的纵向反作用力， 为第i时刻斜率法估算的路面附着系数；
[0106] 43)为了使得代价函数最小，设计基于斜率法的递归最小二乘法为：
[0107]
[0108]
[0109]
[0110] 式中， 表示第k时刻斜率法估算的路面附着系数； 表示第 k‑1时刻‑1
斜率法估算的路面附着系数；κ (k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；Ls(k)表示第k时刻基于斜率法的递推最小二成法增益矩阵；Ps(k)表示第k时刻基于斜率法的最小二乘法更新矩阵；Ps(k‑1)表示第k‑1时刻基于斜率法的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵；I表示单位矩阵；
[0111] 44)通过式(35)‑(37)迭代，估算获得斜率法估算的路面附着系数
[0112] 最终，修正后的路面附着系数估计结果表示为：
[0113]
[0114] 式中， 表示修正后的路面附着系数估计结果；κ表示轮胎滑移率； 表示基于斜率法估计的路面附着系数；表示基于刷子模型估计的路面附着系数。
[0115] 本发明的有益效果为：
[0116] 1)本发明搭建的车辆动力学模型完整的考虑了车辆行驶过程的非线性动力学问题，充分表征了智能电动汽车在结构化道路行驶中所表现的关键行为特征；
[0117] 2)本发明基于卡尔曼滤波算法设计的轮胎力观测器有效地降低了轮胎力观测值的不确定性，为路面附着系数估计提供精确、稳定的轮胎力信息；
[0118] 3)本发明基于最小二乘法和刷子模型能够精确的估计在车辆滑移率较高时的路面附着系数；
[0119] 4)本发明基于最小二乘法和斜率法能够在车辆滑移率较低时对估计的路面附着系数结果进行合理的修正。
附图说明
[0120] 为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
[0121] 图1为本发明的算法构图；
[0122] 图2为三自由度车辆动力学模型示意图；
[0123] 图3为车轮受力示意图；
[0124] 图4为低附着路面轮胎力观测器结果图；
[0125] 图5为低附着路面下路面附着系数估计结果图；
[0126] 图6为高附着路面轮胎力观测器结果图；
[0127] 图7为高附着路面下路面附着系数估计结果图。
具体实施方式
[0128] 下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部结构。
[0129] 参阅图1，一种基于卡尔曼滤波和最小二乘法的车辆路面附着系数估计方法，包括以下步骤：
[0130] 步骤一、参阅图2，为了便于分析，忽略车辆俯仰、侧倾特性，假设车辆左右轮受力对称，基于车辆动力学模型对车辆的状态参数进行求解；
[0131] 所述状态参数包括车速、制动压力、前轮转角、轮胎滑移率、侧偏角和车轮垂向力；所述车速、制动压力、前轮转角由传感器获得；所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算出。
[0132] 所述轮胎滑移率、侧偏角和车轮垂向力由车辆动力学模型计算的具体方法如下：
[0133] 11)建立包含车辆纵向、侧向和横摆运动的三自由度模型，得车辆动力学模型方程为：
[0134] max＝[cosδcosδ1 1 ‑sinδ0]Ftire (1)
[0135] may＝[sinδsinδ0 0 cosδ1]Ftire (2)
[0136]
[0137] 式中，Ftire＝[Fxfl Fxfr Fxrl Fxrr Fyf Fyr]T表示轮胎力集合；Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力；Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyrl表示后左轮侧向力；Fyrr表示后右轮侧向力；m表示整车质量，ax表示车辆纵向加速度；ay表示车辆侧向加速度；δ表示车辆前轮转角；Iz表示车辆绕z轴的转动惯量；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；tw表示轮距；表示车辆质心处横摆角加速度；
[0138] 12)车辆在行驶过程中，车轮在地面上行驶时的受力情况如图3所示。根据动力学平衡关系得到车轮的动力学方程为：
[0139]
[0140] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Iwi表示轮胎旋转惯性质量； 表示车轮转动角加速度；Tdi表示车轮的驱动转矩；Tbi表示车轮的驱动转矩与制动转矩；Fxi表示车轮受到地面的纵向反作用力；r表示车辆滚动半径；TFi＝(a+bvi)Fzir表示车轮受到的滚动阻力矩；a和b表示滚动阻力系数参数值；vi表示车轮中心处的纵向速度；Fzi表示车轮受到地面的垂向反作用力；
[0141] 13)智能电动汽车的四轮驱动力矩可以表示为：
[0142]
[0143] 式中，Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Te表示整车动力电机输出转矩；i0表示传动系的传动比；ηt表示传动系的传动效率；
[0144] 14)车辆制动过程中由制动力形成的制动力矩表示为：
[0145] Tbi＝kbiPwi (6)
[0146] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Tbi表示车轮的制动转矩；kbi表示制动效能因数；Pwi表示轮缸制动压力；
[0147] 15)车轮受到地面的垂向反作用力与车辆加速和减速形成的惯性力存在关系，表示为：
[0148]
[0149] 式中，Fzfl表示前左车轮受到地面的垂向反作用力；Fzfr表示前右车轮受到地面的垂向反作用力；Fzrl表示后左车轮受到地面的垂向反作用力；Fzrr表示后右车轮受到地面的垂向反作用力；m表示整车质量；ax表示车辆纵向加速度； ay表示车辆侧向加速度；lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；hg表示质心高度；tw表示轮距；
[0150] 16)选用前左轮纵向力，前右轮纵向力，后左轮纵向力，后右轮纵向力，前轮总侧向力，后轮总侧向力，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统状态变量，即选用车辆纵向加
速度，车辆侧向加速度，前左车轮转动角速度，前右车轮转动角速度，后左车轮转动角速度，后右车轮转动角速度和车辆质心处横摆角速度为系统输出，
那么得到车辆动力学空间方程为：
[0151]
[0152] 式中，X表示系统状态变量； 表示系统状态变量的微分；u表示系统输入；Φ表示系统状态变量矩阵；B表示系统输入矩阵；Y表示系统输出；H表示系统输出矩阵，且有：
[0153]
[0154]
[0155] Hx＝[cosδcosδ1 1 ‑sinδ0],
[0156] Hy＝[sinδsinδ0 0 cosδ1],
[0157]
[0158] 式中，Φ表示系统状态变量矩阵；B表示系统输入矩阵；H表示系统输出矩阵；r表示车辆滚动半径；I表示单位矩阵；Tdfl表示前左车轮的驱动转矩；Tdfr表示前右车轮的驱动转矩；Tdrl表示后左车轮的驱动转矩；Tdrr表示后右车轮的驱动转矩；Tbfl表示前左车轮的制动转矩；Tbfr表示前右车轮的制动转矩；Tbrl表示后左车轮的制动转矩；Tbrr表示后右车轮的制动转矩；Iwfl表示前左车轮的轮胎旋转惯性质量；Iwfr表示前右车轮的轮胎旋转惯性质量；Iwrl表示后左车轮的轮胎旋转惯性质量；Iwrr表示后右车轮的轮胎旋转惯性质量；Iz表示车辆绕Z轴的转动惯量；Hx表示纵向力矩阵；Hy表示侧向力矩阵； 表示横摆角矩阵；lf表示车辆质心至前轴的距离；lr表示车辆质心至后轴的距离；δ表示车辆前轮转角；tw表示轮距；m表示整车质量；TFfl表示前左车轮受到的滚动阻力矩；TFfr表示前右车轮受到的滚动阻力矩；
TFrl表示后左车轮受到的滚动阻力矩； TFrr表示后右车轮受到的滚动阻力矩；
[0159] 17)车辆轮胎的纵向力和侧向力与路面附着系数紧密相关。采用刷子模型表征车辆轮胎的纵向力和车辆轮胎的侧向力之间的联系，表示为：
[0160]
[0161]
[0162]
[0163]
[0164] 式中，i＝fl,fr,rl,rr分别表示对应的前左车轮，前右车轮，后左车轮和后右车轮；Fxi表示车轮受到地面的纵向反作用力；Fyi表示车轮受到地面的侧向反作用力；Cxi表示轮胎的纵向刚度；Cαi表示轮胎的侧向刚度；κi表示轮胎的滑移率；αi表示轮胎侧偏角；fi表示刷子模型刚度系数；Fti表示刷子模型表征力；μ表示路面附着系数；Fzi表示车轮受到地面的垂向反作用力；
[0165] 18)根据车轮旋转角速度、前轮转角、车速以及车辆结构参数，计算各个轮胎的滑移率：
[0166]
[0167]
[0168]
[0169]
[0170] 式中，κfl表示前左轮滑移率；κfr表示前右轮滑移率；κrl表示后左轮滑移率；κrr表示后右轮滑移率；ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度；r表示车辆滚动半径；vx表示车辆纵向速度；vy表示车辆侧向速度；tw表示轮距；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；lf表示车辆质心到前轴距离；表示车辆质心处横摆角速度；
[0171] 19)智能电动汽车在行驶过程中，前轮转角较小，对应轮胎侧偏角为：
[0172]
[0173] 式中，αfl表示前左轮轮胎侧偏角；αfr表示前右轮轮胎侧偏角；αrl表示后左轮轮胎侧偏角；αrr表示后右轮轮胎侧偏角；表示车辆质心处横摆角速度； lf表示车辆质心到前轴距离；lr表示车辆质心到后轴距离；δfl表示车辆前左轮转角；δfr表示车辆前右轮转角；β表示车辆质心侧偏角；vx表示车辆纵向速度；
[0174] 步骤二、基于卡尔曼滤波估计出车辆行驶过程中轮胎受到的轮胎力；该算法是基于已知的观测值，通过最小化估计值的均方差实现对未知值的估计。
[0175] 具体方法如下：
[0176] 21)将连续的车辆动力学空间方程(8)转换为离散系统状态空间方程，即：
[0177]
[0178] 式中，X(k)表示第k时刻离散系统的状态变量；u(k)表示第k时刻离散系统的输入；X(k+1)表示第k+1时刻离散系统的状态变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；Y(k)表示第k 时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵；w(k)表示第k时刻离散系统的激励噪声；v(k)表示第k时刻离散系统的观测噪声；
[0179] 22)假设离散系统的激励噪声w(k)和观测噪声v(k)是均值为零、方差各为 Q和R的不相关白噪声，离散系统的初始状态X(0)不相关于激励噪声w(k)和观测噪声v(k)；那么卡尔曼滤波的一步预测为：
[0180]
[0181] 式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值； 表示第k时刻的离散系统变量；Φk表示第k时刻离散系统的状态变量矩阵；Bk表示第k时刻离散系统的输入矩阵；u(k) 表示第k时刻离散系统的输入；
[0182] 23)离散系统状态更新矩阵为：
[0183]
[0184] 式中， 表示根据第k时刻离散系统的状态变量对k+1时刻离散系统的状态变量进行预测的值；K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵； Y(k+1)表示第k+1时刻离散系统的输出；Hk表示第k时刻离散系统的输出矩阵； 表示经校正后第k+1
时刻卡尔曼滤波的估计值；
[0185] 24)卡尔曼滤波增益矩阵表示为：
[0186]
[0187] 式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Hk表示第k 时刻离散系统的输出矩阵；R表示观测噪声v(k)的方差；
[0188] 25)一步预测协方差矩阵为：
[0189]
[0190] 式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第 k时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；
[0191] 26)在一步预测协方差矩阵中引入遗忘因子进行修正，即：
[0192]
[0193] 式中，P(k+1|k)表示根据第k时刻的误差协方差与状态向量对k+1时刻的预测协方差矩阵；Φk表示第k时刻离散系统的状态变量矩阵；P(k|k)表示第k 时刻离散系统的协方差矩阵；Q表示激励噪声w(k)的方差；s表示遗忘因子；
[0194] 协方差矩阵更新方程为：
[0195] P(k+1|k+1)＝(I‑K(k+1)Hk)P(k+1|k) (24)
[0196] 式中，K(k+1)表示第k+1时刻卡尔曼滤波增益矩阵；Hk表示第k时刻离散系统的输出矩阵；P(k+1|k)表示根据第k时刻的误差协方差与状态向量对 k+1时刻的预测协方差矩阵；P(k+1|k+1)为经校正后的协方差矩阵；I表示单位矩阵；
[0197] 27)通过式 (19)‑ (2 4) 迭代 ，估算获得 离散系统的状态变 量中的轮胎纵向力和
侧向力；
[0198] 其中，Fxfl表示前左轮纵向力；Fxfr表示前右轮纵向力；Fxrl表示后左轮纵向力；Fxrr表示后右轮纵向力；Fyf＝(Fyfl+Fyfr)表示前轮总侧向力； Fyr＝(Fyrl+Fyrr)表示后轮总侧向力；Fyfl表示前左轮侧向力；Fyfr表示前右轮侧向力；Fyfl表示后左轮侧向力；Fyrr表示后右轮侧向力；ωfl表示前左轮转动角速度；ωfr表示前右轮转动角速度；ωrl表示后左轮转动角速度；ωrr表示后右轮转动角速度，表示车辆质心处横摆角速度。
[0199] 步骤三、采用最小二乘法对路面附着系数进行了估计；
[0200] 具体方法如下：
[0201] 31)对刷子模型进行改写，即：
[0202] y(k)＝f(k,θ(k)) (25)
[0203] 式中，y(k)＝[Fx,Fy]表示第k时刻的轮胎力；Fx表示车轮受到地面的纵向反作用力；Fy表示车轮受到地面的侧向反作用力；f(k,θ(k))表示第k时刻的刷子轮胎模型；θ(k)＝T[Cx,Cα,μ]表示第k时刻刷子轮胎模型的状态向量；Cx表示轮胎的纵向刚度；Cα表示轮胎的侧向刚度；μ表示路面附着系数；
[0204] 32)对y(k)进行一阶泰勒展开，并忽略高阶项，得：
[0205] y(k)≈F(k)(θ(k)‑θ(k‑1))+f(θ(k‑1),k‑1) (26)
[0206] 式中，y(k)表示第k时刻的轮胎力；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；θ(k)为第k时刻刷子轮胎模型的状态向量；f(θ(k‑1),k‑1)表示第k‑1时刻的刷子轮胎模型；
F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置，且有：
[0207]
[0208] 式中，F(k)表示y(k)对θ(k)的雅可比矩阵的转置；θ(k)为第k时刻刷子轮胎模型的状态向量；f(k,θ)表示刷子轮胎模型；θ表示刷子轮胎模型的状态向量；θ(k‑1)为第k‑1时刻刷子轮胎模型的状态向量；
[0209] 令z(k)＝y(k)‑f(θ(k‑1),k‑1)+F(k)θ(k‑1)，那么：
[0210] z(k)≈F(k)θ(k) (28)
[0211] 式中，F(k)表示第k时刻y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k 时刻的轮胎力；θ(k)第k时刻刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；
[0212] 33)通过最小二乘法对路面附着系数进行估计，设定最小二乘法代价函数为：
[0213]
[0214] 式中， 表示基于刷子轮胎模型的递归最小二乘法代价函数；Λ表示最小二乘法的遗忘因子矩阵；z(i)表示第i时刻可测量轮胎力矩阵；F(i)表示第i时刻y(i)对θ(i)的雅可比矩阵的转置；θ(i)为第i时刻刷子轮胎模型的状态向量； 表示第i时刻估计的刷子轮胎模型的状态向量；
[0215] 34)为了使得代价函数最小，基于刷子轮胎模型的递归最小二乘法为：
[0216]
[0217] L(k)＝P(k‑1)FT(k)(I+F(k)P(k‑1)FT(k))‑1 (31)
[0218] P(k)＝Λ‑1(I‑L(k)F(k))P(k‑1)Λ‑1 (32)
[0219] 式中， 表示第k时刻估计的刷子轮胎模型的状态向量； 表示第k‑1时刻估计的刷子轮胎模型的状态向量；z(k)为第k时刻可测量轮胎力矩阵；F(k‑1)表示第k‑1时刻y(k‑1)对θ(k‑1)的雅可比矩阵的转置；y(k‑1) 表示第k‑1时刻的轮胎力，θ(k‑1)第k‑1时刻刷子轮胎模型的状态向量，L(k) 表示第k时刻基于刷子模型的递推最小二成法增益矩阵；F(k)表示第k时刻 y(k)对θ(k)的雅可比矩阵的转置；y(k)表示第k时刻的轮胎力；P(k)表示第k 时刻基于刷子模型的最小二乘法更新矩阵；P(k‑1)表示第k‑1时刻基于刷子模型的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵，I表示单位矩阵； 表示××；
[0220] 35)通过(30)‑(32)迭代，估算获得估计的刷子轮胎模型的状态向量中的路面附着系数；其中， 表示基于刷子模型估计的轮胎的纵向刚度，
表示基于刷子模型估计的轮胎的侧向刚度，表示基于刷子模型估计的路面附着系数。
[0221] 步骤四、通过斜率法对步骤三中低滑移率下的路面附着系数估计结果进行修正。
[0222] 具体方法如下：
[0223] 41)在滑移率小时，附着系数μ与滑移率κ近似为正比例关系，即：
[0224]
[0225] 式中，κ‑1(k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；μs(k)表示第k时刻用于斜率法真实的路面附着系数；
[0226] 42)设定基于斜率法的最小二乘法代价函数为：
[0227]
[0228] 式中， 表示基于斜率法的递归最小二乘法代价函数；Λ表示最小二乘法的‑1遗忘因子矩阵；κ (i)表示第i时刻轮胎的滑移率的倒数；Fz(i)表示第i时刻车轮受到地面的垂向反作用力，Fx(i)表示第i时刻车轮受到地面的纵向反作用力， 为第i时刻斜率法估算的路面附着系数；
[0229] 43)为了使得代价函数最小，设计基于斜率法的递归最小二乘法为：
[0230]
[0231]
[0232]
[0233] 式中， 表示第k时刻斜率法估算的路面附着系数； 表示第 k‑1时刻‑1
斜率法估算的路面附着系数；κ (k)表示第k时刻轮胎的滑移率的倒数；Fz(k)表示第k时刻车轮受到地面的垂向反作用力；Fx(k)表示第k时刻车轮受到地面的纵向反作用力；Ls(k)表示第k时刻基于斜率法的递推最小二成法增益矩阵；Ps(k)表示第k时刻基于斜率法的最小二乘法更新矩阵；Ps(k‑1)表示第k‑1时刻基于斜率法的最小二乘法更新矩阵；Λ表示最小二乘法的遗忘因子矩阵；I表示单位矩阵；
[0234] 44)通过式(35)‑(37)迭代，估算获得斜率法估算的路面附着系数
[0235] 最终，修正后的路面附着系数估计结果表示为：
[0236]
[0237] 式中， 表示修正后的路面附着系数估计结果；κ表示轮胎滑移率； 表示基于斜率法估计的路面附着系数；表示基于刷子模型估计的路面附着系数。
[0238] 实施例
[0239] 我们基于MATLAB/Simulink与车辆动力学软件CarSim搭建了联合仿真平台用于对本申请的轨迹跟踪控制器进行测试。
[0240] 第一组工况为低附着路面下的制动工况，设置0s时车辆制动系统开始介入，且制动压力为2Mpa，路面附着系数为0.3，车辆初始车速为100km/h。实验结果如图4、图5所示。
[0241] 从图中可以看出，所建立的轮胎力观测器在低附着路面上轮胎力和路面附着系数估计精度较高。路面附着系数估计算法能够较快的收敛到估计的路面附着系数，约为0.2816左右。与道路实际路面附着系数相比，估计的相对误差为6.13％，满足实际车辆运动控制使用需求。
[0242] 第二组工况为高附着路面下的制动工况，设置路面的附着系数为0.7，并且设置车辆的制动压力为2MPa，制动起始时间为0s，车辆初始车速为100km/h。实验结果如图6、图7所示。
[0243] 从图中可以看出，测器收敛后轮胎纵向力的估计误差最大值为202.76N，此时估计的相对误差6.31％。路面附着系数估计值为0.7152，估计的相对误差2.17％，同时RLS算法能在0.4s内收敛，说明高附着路面下识别算法具有不错的响应速度。
