CN104691551B - 一种基于电机与车轮耦合特性的路面附着系数估计方法 - Google Patents

Abstract

本发明公开一种基于电机与车轮耦合特性的路面附着系数估计方法，适用于电动车辆行驶过程中路面附着系数的实时监测。它建立轮胎纵向刚度与电动轮共振频率之间的关系，在利用电机转矩获取的共振频率基础上求解出轮胎纵向刚度，然后再利用轮胎纵向刚度与路面附着系数的关系，实现了路面附着系数估计。本发明仅采用电机电流与轮速信号，不需要车速与轮胎纵向力信息，不需要计算轮胎纵向滑移率，使得该方法应用方便；利用频域信息进行估计，使得该方法具有对轮速噪声与误差不敏感的特性，也说明了该方法的准确性。

Description

一种基于电机与车轮耦合特性的路面附着系数估计方法

技术领域

本发明涉及车辆的路面附着系数估计，特别是关于一种分布式驱动的电动车辆的路面附着系数估计方法，是基于电机与车轮耦合特性的路面附着系数估计。

背景技术

分布式电驱动车辆是将驱动电机分别安装在各车轮内或各车轮附近，一台电机独立驱动一个车轮，具有响应速度快、传动链短、传动高效、结构紧凑等优点，其充分利用了电机转矩精确观测和快速可控的特点，是电动汽车领域的一个重要发展方向。但目前的分布式电驱动车辆仍存在很多需要改进的问题，如低速时，电机转矩波动造成的冲击，对于该现象还没有合理的解释。

路面附着系数：是指轮胎与地面间作用的纵向力、侧向力的合力与垂向力之比的最大值。精确估计路面附着系数是研究电机转矩波动造成冲击影响的可靠前提。目前国内外对于路面峰值附着系数实时估算方法已经进行了大量研究。这些方法可以分为基于原因的方法和基于效果的方法两类。前者是利用超声波传感器等来检测路面状况来估算路面附着系数，该种方法需要外加昂贵的传感器,并且对于环境的依赖程度较高。后者方法则是直接利用车辆与轮胎的动力学特性来估计路面附着系数，例如用μ-s曲线斜率(附着系数与滑移率曲线)估算路面附着系数的方法。该类方法由于需要准确的纵向力和滑移率估计值，所以对轮速噪声和稳态误差的要求比较高，也就是需要车轮发生较大滑转时才能较好的工作，而且目前该类方法均采用稳态轮胎模型，不适用于瞬态工况，特别是在分布式电驱动车辆这种结构下，高频振动源较多，使得纵向力不能通过稳态轮胎模型估计。

发明内容

为了解决现有用μ-s曲线斜率估算路面附着系数的方法不适用于小滑移率工况，且对轮速噪声和稳态误差灵敏度高的缺点，本发明提出一种新型的基于电机与车轮耦合特性的路面附着系数估计方法，适用于电驱动车辆正常行驶(匀速或者小幅加/减速)过程中路面附着系数的实时监测估计，可在不需要转矩传感器和车速传感器等情况下，仅依靠轮速信号处理实现对每个车轮所在路面情况的辨识。

为实现上述目的，本发明采取以下技术方案：一种基于电机与车轮耦合特性的路面附着系数估计方法，其特征在于，包括以下步骤：

1)首先建立轮胎纵向刚度与车轮共振频率之间的关系，为：

$f_0\≈\frac{R}{2\mathrm{\π}}\sqrt{\frac{k_s}{{\mathrm{Ir}}_x}}---\left(20\right)$

式中f₀为共振频率，R为车轮滚动半径，k_s为车轮纵向刚度，I为车轮转动惯量，r_x为轮胎的纵向松弛长度；

2)然后在整车控制器获取实时的轮速信号ω和电机的电流信号i_q基础上，利用MATLAB中的nonlinear ARX model模块，输出二阶系统模型传递函数的系数a₁，a₂，a₃，然后利用公式(22)找到两个解λ_i(i＝1,2)：

G(λ)＝a₁λ²+a₂λ+a₃ (22)

式中λ是拉布拉斯算子，a₁，a₂，a₃为系数；

3)再按照公式(23)、(24)、(25)计算共振频率f₀：

$f_0=\frac{\sqrt{d_i^2-c_i^2}}{2\mathrm{\π}}---\left(23\right)$

式中 $c_i=\frac{\ln (\mathrm{Re}{\left({\mathrm{\λ}}_i\right)}^2+\mathrm{Im}{\left({\mathrm{\λ}}_i\right)}^2)}{2\mathrm{\ΔT}}---\left(24\right)$

$d_i=-\frac{1}{2\mathrm{\ΔT}}\frac{\mathrm{Im}\left({\mathrm{\λ}}_i\right)}{\mathrm{Re}\left({\mathrm{\λ}}_i\right)}---\left(25\right)$

上述ΔT为采样时间，Re、Im分别表示数学计算中的实部和虚部；

4)将公式(23)取得的f₀带入到公式(20)中，在公式(20)的基础上估计轮胎纵向刚度k_s；

5)在任意时刻，将得到的轮胎纵向刚度输入到路面附着系数估计模块中，计算得到路面附着系数：

${\mathrm{\μ}}_{\max }\left(t\right)=a_{{\mathrm{\μ}}_{\max }}{k}_s\left(t\right)+b_{{\mathrm{\μ}}_{\max }}---\left(25\right)$

式中的系数根据实验数据确定，采用数据拟合的方法得到。

在建立轮胎纵向刚度与车轮共振频率之间的关系时，是按照如下方法实现的：

1)建立单轮动力学模型：

$I\stackrel{\·}{\mathrm{\ω}}=T_d-{F_x}^{D}R---\left(1\right)$

式中I为车轮转动惯量，ω为车轮转速，是ω关于时间的导数，T_d为电机输出转矩，F_x ^D是车轮瞬态纵向力，R为车轮滚动半径；

2)建立稳态轮胎模型：

$F_x^S=k_{s}s+F_{x0}$

式中是车轮稳态纵向力，k_s为车轮纵向刚度，s为滑移率，不同情况取不同值，F_x0为滑移率为零时的纵向力，v为车辆纵向速度，ω为车轮转速，R为车轮滚动半径；

3)建立瞬态轮胎模型：

$\begin{array}{c}\mathrm{\τ}{{\stackrel{\·}{F}}_x}^D+{F_x}^D={F_x}^S\\ \mathrm{\τ}=\frac{r_x}{v}\end{array}---\left(3\right)$

式中F_x ^D为车轮瞬态纵向力，是F_x ^D关于时间的导数，τ为时间常数，为车轮稳态纵向力，r_x为轮胎的纵向松弛长度；

4)建立简化电机模型：

T_d＝Ki_q (4)

式中i_q为电机电流，K为比例常数，通过实验测得或电机厂商提供，T_d为电机输出转矩；

5)在上述几种车辆模型和简化电机模型的基础上，假设电机转矩由两部分组成，分为恒定部分与高频部分，如下所示：

T_d＝T₀+T₁sin(2πf·t) (5)

式中T₀表示恒定转矩，认为是一个相对恒定的值；而T₁sin(2πf·t)表示高频转矩，其中T₁为转矩振幅，f为高频转矩的频率，t是表示某一时刻；

6)结合公式(1)，将公式(5)带入到公式(1)中，单轮动力学模型可表示为：

$I\stackrel{\·}{\mathrm{\ω}}=T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{D}R---\left(6\right)$

对公式(6)两端同时对时间求导，得到公式(7)：

$I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\πf}\·T_1\cos (2\mathrm{\πf}\·t)-{\stackrel{\·}{F}}_x^{D}R---\left(7\right)$

再将公式(7)乘以时间常数τ，然后与公式(6)求和得公式(8)：

$\begin{array}{c}\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{D}R\\ +2\mathrm{\πf\τ}\·T_1\cos (2\mathrm{\πf}\·t)-\mathrm{\τ}{\stackrel{\·}{F}}_x^{D}R\\ =T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{S}R\\ +2\mathrm{\πf\τ}\·T_1\cos (2\mathrm{\πf}\·t)\end{array}---\left(8\right)$

结合公式(2)、(3)和三角函数公式，公式(8)可简化为公式(9)：

$\begin{array}{c}\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+a{T}_1\sin (2\mathrm{\πf}\·t+\mathrm{\φ})-(k_{s}s+F_{x0})R\\ a=\sqrt{1+{\left(2\mathrm{\πf}\right)}^2},\mathrm{\φ}=\arctan \left(2\mathrm{\πf}\right)\end{array}---\left(9\right)$

式中aT₁为合并后高频信号的振幅，φ合并后高频信号的初始相位；

对于驱动工况，考虑其滑移率定义，公式(9)可以表示为公式(10)：

$\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+a{T}_1\sin (2\mathrm{\πf}\·t+\mathrm{\φ})-k_{s}R\frac{\mathrm{R\ω}-v}{\mathrm{R\ω}}-F_{x0}R---\left(10\right)$

进一步对公式(10)两端求导可得公式(11)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=a2\mathrm{\πf}\·T_1\cos (2\mathrm{\πf}\·t+\mathrm{\φ})+k_{s}R\frac{\stackrel{\·}{v}\·\mathrm{R\ω}-R\stackrel{\·}{\mathrm{\ω}}\·v}{{\left(\mathrm{R\ω}\right)}^2}---\left(11\right)$

假设Rω≈v，由于车辆的惯性远大于车轮的惯性，所以车辆的加速度相比车轮的角加速度可以忽略，公式(11)可进一步简化为公式(12)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\πfa}\·T_1\cos -(2\mathrm{\πf}\·t+\mathrm{\φ})-\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}---\left(12\right)$

令：

T₂＝T₁cos(2πf·t+φ) (13)

公式(12)可表示为公式(14)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}+\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}=2\mathrm{\πfa}\·T_2---\left(14\right)$

对公式(14)两端作拉普拉斯变换，得到公式(15)：

$\mathrm{I\τ}{\mathrm{\λ}}^3\mathrm{\ω}\left(\mathrm{\λ}\right)+I{\mathrm{\λ}}^2\mathrm{\ω}\left(\mathrm{\λ}\right)+\frac{k_s{R}^2}{v}\mathrm{\λ\ω}\left(\mathrm{\λ}\right)=2\mathrm{\πfa}{T}_2\left(\mathrm{\λ}\right)---\left(15\right)$

式中λ为拉普拉斯算子；

这样就得到了电机转矩到轮速的传递函数，公式(16)：

$\frac{T_2\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{\mathrm{I\τ}{\mathrm{\λ}}^3+I{\mathrm{\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\πf}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}---\left(16\right)$

结合公式(4)和(16)，可以进一步得到可以得到电机电流到轮速的传递函数，公式(17)：

$\frac{i_q\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{\mathrm{I\τ}{\mathrm{\λ}}^3+I{\mathrm{\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\πfK}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}---\left(17\right)$

令λ＝j2πf，其中j表示虚部，合并同类项，求模即可得到电机电流到轮速的幅频函数，公式(18)：

$\begin{array}{c}\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}=\left|\frac{i_q\left(j2\mathrm{\πf}\right)}{\mathrm{\ω}\left(j2\mathrm{\πf}\right)}\right|\\ =\frac{|-I{\left(2\mathrm{\πf}\right)}^2+j(\frac{k_s{R}^2}{v}2\mathrm{\πf}-\mathrm{I\τ}{\left(2\mathrm{\πf}\right)}^3)|}{2\mathrm{\πfK}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\\ =\frac{|-I\left(2\mathrm{\πf}\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_s{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\end{array}---\left(18\right)$

作如下近似：当电机车轮系统发生共振时，即有公式(19)：

$\begin{array}{c}\min \left[\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}\right]=\min \left[\frac{|-I\left(2\mathrm{\πf}\right)+\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\right]\\ \≈\left[\frac{|-I\left(2\mathrm{\πf}\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\·2\mathrm{\πf}}\right]\\ =\min \left[\right|\frac{-I}{K}+j\left(\frac{k_s{R}^2/2\mathrm{\πf}-{\mathrm{Ir}}_x\left(2\mathrm{\πf}\right)}{\mathrm{vK}}\right)\left|\right]\end{array}---\left(19\right)$

应用求极值的方法，得到其最小值对应的频率，即共振频率，见公式(20)。

由于轮胎纵向刚度是随时间变化的，所以基于递归最小二乘法进行轮胎纵向刚度估计，递归最小二乘法的轮胎纵向刚度估计可以描述为以下方程：

x(t)＝θ^T(t)k_s(t)+e(t) (26)

其中k_s(t)是随时间变化的轮胎纵向刚度；x(t)是共振频率f₀的平方；e(t)是误差函数，θ^T(t)是矩阵θ(t)的转置，矩阵是回归系数函数，θ(t)为式中R为车轮滚动半径，I为车轮转动惯量，r_x为轮胎松弛长度；

对式(26)求解k_s(t)按如下步骤进行：

步骤一：获取系统输出的x(t)，并确定回归系数函数θ^T(t)；

步骤二：根据t时刻的系统真实输出值和上一时刻t-1预测的本时刻输出值来计算t时刻e(t)：

e(t)＝x(t)-θ^T(t)k_s(t-1) (27)

步骤三：利用下式求t时刻解增益向量K(t):

$K\left(t\right)=\frac{P(t-1)\mathrm{\θ}\left(t\right)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}---\left(28\right)$

式中P(t-1)为上一时刻t-1预测的本时刻的协方差矩阵，参见式(29)，χ是遗忘因子,本例中取值在(0.9，1)，

利用下式计算下一时刻t+1的协方差矩阵P(t):

$P\left(t\right)=\frac{1}{\mathrm{\χ}}[P(t-1)-\frac{P(t-1)\mathrm{\θ}\left(t\right){\mathrm{\θ}}^T\left(t\right)P(t-1)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}]---\left(29\right)$

步骤四：计算轮胎纵向刚度：

k_s(t)＝k_s(t-1)+K(t)e(t) (30)。

对于式(26)的求解可按如下步骤进行：

步骤一：获取系统输出的x(t)，即上一模块估计共振频率的平方，并确定回归系数函数θ^T(t)。

步骤二：根据t时刻的系统真实输出值和上一时刻t-1预测的本时刻输出值来计算t时刻e(t)：

e(t)＝x(t)-θ^T(t)k_s(t-1) (27)

步骤三：利用下式求t时刻解增益向量K(t):

$K\left(t\right)=\frac{P(t-1)\mathrm{\θ}\left(t\right)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}---\left(28\right)$

式中P(t-1)为上一时刻t-1预测的本时刻的协方差矩阵，参见式(29)。χ是遗忘因子,本例中取值在(0.9，1)。

利用下式计算下一时刻t+1的协方差矩阵P(t):

$P\left(t\right)=\frac{1}{\mathrm{\χ}}[P(t-1)-\frac{P(t-1)\mathrm{\θ}\left(t\right){\mathrm{\θ}}^T\left(t\right)P(t-1)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}]---\left(29\right)$

步骤四：计算轮胎纵向刚度：

k_s(t)＝k_s(t-1)+K(t)e(t) (30)。

本发明所展现出来的优点是：

1、本方法推导得到的电机车轮共振频率的公式能够表征汽车正常行驶工况(匀速、小加(减)速)下轮胎和路面特征参数对电机车轮系统动态响应的影响。

2、本方法在电机车轮共振频率的公式的基础上，仅采用电机电流与轮速信号，不需要车速与轮胎纵向力信息，不需要添加额外的传感器，不需要计算轮胎纵向滑移率，即可以估计出轮胎纵向刚度，进而得到路面附着系数，应用方便。

3、本方法利用频域信息进行估计，使得该方法具有对轮速噪声与误差不敏感的特性，也说明了该方法的准确性。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分的从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

附图仅用于示出具体实施例的目的，并不是对本发明的限制。

图1是本发明方法的实施流程图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述，其中，附图用于与本发明的实施例一起用于阐释本发明，但本领域的技术人员应该知道，以下实施例并不是对本发明技术方案作的唯一限定，凡是在本发明技术方案精神实质下所做的任何等同变换或改动，均应视为属于本发明的保护范围。

本发明是在大量的理论与实践相结合的基础上，首先摸索建立轮胎纵向刚度与电动轮共振频率之间的关系，然后通过此两者的关系，依靠一套植有计算算法的估计系统来实现路面附着系数的估计。该系统包括：一设置在车轮处的轮速传感器，设置在整车控制器中的一轮边电机控制器、一基于nonlinear ARX model模型(非线性ARX模型)的电机车轮共振频率估计模块、一附着系数估计模块、一递归最小二乘法的轮胎纵向刚度估计模块，这些模块都附有算法程序，集成在整车控制器的控制系统中。

因为常规车辆车轮的转矩是不可控的，所以不能够通过控制车轮转矩使其发生共振现象，而电动车辆在不同频率下的转矩是可控的，所以可以通过控制车轮转矩使车轮发生共振，共振特性存在于电动轮中，所以我们可以用这个方法对电动车辆的路面附着系数进行估计。

一，本发明首先建立轮胎纵向刚度与电动轮共振频率之间的关系，其采取了如下方法：

1、首先建立如下几个车辆模型：

1)建立单轮动力学模型：

$I\stackrel{\·}{\mathrm{\ω}}=T_d-{F_x}^{D}R---\left(1\right)$

这个计算模型表达式是公知的。

式中I为车轮转动惯量，ω为车轮转速，是ω关于时间的导数，T_d为电机输出转矩，同时也是车轮驱动转矩，F_x ^D是车轮瞬态纵向力，R为车轮滚动半径。这些参数中，I、R是车辆已知的，ω是轮速传感器获得的，T_d是通过电机的电流换算出来的，见下面公式(4)。

2)建立稳态轮胎模型：

$F_x^S=k_{s}s+F_{x0}---\left(2\right)$

上述表达式也是公知的。

式中是车轮稳态纵向力，k_s为车轮纵向刚度，s为滑移率，不同情况取不同值，F_x0为滑移率为零时的纵向力，v为车辆纵向速度，ω为车轮转速，这些也都为已知或可求，但是在本发明中作为中间量不用具体取值。

3)建立瞬态轮胎模型：

$\begin{array}{c}\mathrm{\τ}{{\stackrel{\·}{F}}_x}^D+{F_x}^D={F_x}^S\\ \mathrm{\τ}=\frac{r_x}{v}\end{array}---\left(3\right)$

上述表达式也是公知的。

式中F_x ^D为车轮瞬态纵向力，是F_x ^D关于时间的导数，τ为时间常数，为车轮稳态纵向力，r_x为轮胎的纵向松弛长度，松弛长度可由研究车辆模型的人提供，为已知量。

2、其次建立简化电机模型：

T_d＝Ki_q (4)

这一模型表达式也是已知的。

式中i_q为电机电流，是系统直接获得的，K为比例常数，通过实验测得或电机厂商提供，T_d为电机输出转矩。

3、在上述几种车辆模型和简化电机模型的基础上，我们进一步推展如下：

我们假设电机转矩由两部分组成，分为恒定部分与高频部分，如下所示：

T_d＝T₀+T₁sin(2πf·t) (5)

式中T₀表示恒定转矩，可以认为是在驾驶员踩踏板时，根据驾驶员需求系统所换算出来的一个转矩，它随时间变化较慢，所以可以认为是一个相对恒定的值；而高频转矩我们一般用T₁sin(2πf·t)表示，这一部分是得到驾驶员需求转矩之后，为了实现我们的这个附着系数估计，加入到电机控制器中的一个高频转矩，其中T₁为转矩振幅，f为高频转矩的频率；t是表示某一时刻。

结合公式(1)，将公式(5)带入到公式(1)中，单轮动力学模型可表示为公式(6)：

$I\stackrel{\·}{\mathrm{\ω}}=T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{D}R---\left(6\right)$

为了用轮胎稳态力表达轮胎瞬态力，我们对公式(6)两端同时对时间求导，得到公式(7)：

$I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\πf}\·T_1\cos (2\mathrm{\πf}\·t)-{\stackrel{\·}{F}}_x^{D}R---\left(7\right)$

再将公式(7)乘以时间常数τ，然后与公式(6)求和得公式(8)：

$\begin{array}{c}\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{D}R\\ +2\mathrm{\πf\τ}\·T_1\cos (2\mathrm{\πf}\·t)-\mathrm{\τ}{\stackrel{\·}{F}}_x^{D}R\\ =T_0+T_1\sin (2\mathrm{\πf}\·t)-F_x^{S}R\\ +2\mathrm{\πf\τ}\·T_1\cos (2\mathrm{\πf}\·t)\end{array}---\left(8\right)$

结合公式(2)、(3)和三角函数公式，公式(8)可简化为公式(9)：

$\begin{array}{c}\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+a{T}_1\sin (2\mathrm{\πf}\·t+\mathrm{\φ})-(k_{s}s+F_{x0})R\\ a=\sqrt{1+{\left(2\mathrm{\πf}\right)}^2},\mathrm{\φ}=\arctan \left(2\mathrm{\πf}\right)\end{array}---\left(9\right)$

式中aT₁为合并后高频信号的振幅，φ合并后高频信号的初始相位。

对于驱动工况，考虑其滑移率定义，公式(9)可以表示为公式(10)：

$\mathrm{I\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+a{T}_1\sin (2\mathrm{\πf}\·t+\mathrm{\φ})-k_{s}R\frac{\mathrm{R\ω}-v}{\mathrm{R\ω}}-F_{x0}R---\left(10\right)$

为了消除T₀和F_x0等常数项，进一步对公式(10)两端求导可得公式(11)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=a2\mathrm{\πf}\·T_1\cos (2\mathrm{\πf}\·t+\mathrm{\φ})+k_{s}R\frac{\stackrel{\·}{v}\·\mathrm{R\ω}-R\stackrel{\·}{\mathrm{\ω}}\·v}{{\left(\mathrm{R\ω}\right)}^2}---\left(11\right)$

通过这些操作就找出了高频电机转矩信息与轮速信息的关系了，但式(11)中仍含有速度的导数项且轮速处于分母中，需要进一步变换才能得到高频电机转矩信息与轮速信息在频域下的关系。

因为以前人们不能够解决滑移率较小工况下的附着系数估计，而本发明方法中不需要估计滑移率，所以适用于滑移率极小情况，所以可以假设Rω≈v，而且车辆的惯性远大于车轮的惯性，即车辆的加速度相比车轮的角加速度可以忽略，在这两个假设前提下，公式(11)可进一步简化为公式(12)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\πfa}\·T_1\cos -(2\mathrm{\πf}\·t+\mathrm{\φ})-\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}---\left(12\right)$

经过以上简化，速度的导数项被忽略。

令：

T₂＝T₁cos(2πf·t+φ) (13)

公式(12)可表示为公式(14)：

$\mathrm{I\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}+\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}=2\mathrm{\πfa}\·T_2---\left(14\right)$

对公式(14)两端作拉普拉斯变换，得到公式(15)，拉普拉斯变换是将时域的表达转换到频域的表达方法，转换后能够找到频域下的振幅和相位信息：

$\mathrm{I\τ}{\mathrm{\λ}}^3\mathrm{\ω}\left(\mathrm{\λ}\right)+I{\mathrm{\λ}}^2\mathrm{\ω}\left(\mathrm{\λ}\right)+\frac{k_s{R}^2}{v}\mathrm{\λ\ω}\left(\mathrm{\λ}\right)=2\mathrm{\πfa}{T}_2\left(\mathrm{\λ}\right)---\left(15\right)$

式中λ为拉普拉斯算子。

这样就得到了电机转矩到轮速的传递函数，如公式(16)：

$\frac{T_2\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{\mathrm{I\τ}{\mathrm{\λ}}^3+I{\mathrm{\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\πf}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}---\left(16\right)$

在实际应用中，能直接测量的信号为电机电流信号，而非电机转矩信号。结合公式(4)和(16)，可以进一步得到可以得到电机电流到轮速的传递函数，如公式(17)：

$\frac{i_q\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{\mathrm{I\τ}{\mathrm{\λ}}^3+I{\mathrm{\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\πfK}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}---\left(17\right)$

令λ＝j2πf，其中j表示虚部，合并同类项，求模即可得到电机电流到轮速的幅频函数如公式(18)所示：

$\begin{array}{c}\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}=\left|\frac{i_q\left(j2\mathrm{\πf}\right)}{\mathrm{\ω}\left(j2\mathrm{\πf}\right)}\right|\\ =\frac{|-I{\left(2\mathrm{\πf}\right)}^2+j(\frac{k_s{R}^2}{v}2\mathrm{\πf}-\mathrm{I\τ}{\left(2\mathrm{\πf}\right)}^3)|}{2\mathrm{\πfK}\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\\ =\frac{|-I\left(2\mathrm{\πf}\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_s{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\end{array}---\left(18\right)$

通过经验数据可知共振频率绝大多数情况下大于10Hz，那么可以作如下近似，当电机车轮系统发生共振时，即有公式(19)：

$\begin{array}{c}\min \left[\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}\right]=\min \left[\frac{|-I\left(2\mathrm{\πf}\right)+\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\πf}\right)}^2}}\right]\\ \≈\left[\frac{|-I\left(2\mathrm{\πf}\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\πf}\right)}^2}{v}\right)|}{K\·2\mathrm{\πf}}\right]\\ =\min \left[\right|\frac{-I}{K}+j\left(\frac{k_s{R}^2/2\mathrm{\πf}-{\mathrm{Ir}}_x\left(2\mathrm{\πf}\right)}{\mathrm{vK}}\right)\left|\right]\end{array}---\left(19\right)$

应用求极值的方法，得到其最小值对应的频率，共振频率的表达式为公式(20)：

$f_0\≈\frac{R}{2\mathrm{\π}}\sqrt{\frac{k_s}{{\mathrm{Ir}}_x}}---\left(20\right)$

上述，前期利用车辆的一系列模型公式，建立电机转矩转换到轮速的动力学方程，然后应用拉普拉斯变换，得到轮速到电机转矩到的传递函数，进一步得到电机电流到轮速的传递函数，具体为电机电流到轮速的幅频传递函数，最后得到车轮纵向刚度k_s、车轮转动惯量I与共振频率f₀之间的关系，这些是我们独创的推演过程。

二，在建立了轮胎纵向刚度与电动轮共振频率之间关系式(式(20))的基础上，开展路面附着系数的估计：

1、整车控制器取实时的轮速信号ω和电机的电流信号i_q，将两个信号发送到MATLAB中的nonlinear ARX model模块(一款公知软件，可参照文献T.Hirao,et al.,Resonance frequency estimation of time-series data by subspace method,Proceedings of the ICROS-SICE International Joint Conference 2009,ISBN9784907764333,4913-4916)，该模块输出二阶系统模型传递函数的系数a₁，a₂，a₃；

G(λ)＝a₁λ²+a₂λ+a₃ (21)

上式中λ是拉布拉斯算子，a₁，a₂，a₃为系数。

2、找到公式(21)的两个解λ_i(i＝1,2)，然后按照公式(22)、(23)、(24)计算共振频率f₀，然后带入到公式(20)中，在公式(20)的基础上估计轮胎纵向刚度。

$f_0=\frac{\sqrt{d_i^2-c_i^2}}{2\mathrm{\π}}---\left(22\right)$

式中 $c_i=\frac{\ln (\mathrm{Re}{\left({\mathrm{\λ}}_i\right)}^2+\mathrm{Im}{\left({\mathrm{\λ}}_i\right)}^2)}{2\mathrm{\ΔT}}---\left(23\right)$

$d_i=-\frac{1}{2\mathrm{\ΔT}}\frac{\mathrm{Im}\left({\mathrm{\λ}}_i\right)}{\mathrm{Re}\left({\mathrm{\λ}}_i\right)}---\left(24\right)$

上述ΔT为采样时间，Re、Im分别表示数学计算中的实部和虚部。

在实际应用中，我们不可能随机获取大量的轮速信号ω和电机的电流信号i_q来进行计算，因此依据前面公式(5)的假设，在车辆实际运行过程中，在驾驶员需求的恒定转矩信号T₀上，我们附加一组有限的高频转矩信号T₃(这组信号为共振频率附近的信号，包含共振频率)，使轮速信号和电流信号在这个频率段下的能量加强，这样会使得估计效果更好，效率更高，所取个数不受限制。

3、在得到轮胎纵向刚度基础上，再进一步估计附着系数：

${\mathrm{\μ}}_{\max }\left(t\right)=a_{{\mathrm{\μ}}_{\max }}{k}_s\left(t\right)+b_{{\mathrm{\μ}}_{\max }}---\left(25\right)$

式中的系数需要根据实验数据确定，采用数据拟合的方法得到。

由于轮胎纵向刚度是时间t的函数，所以可基于递归最小二乘法进行轮胎纵向刚度估计，递归最小二乘法的轮胎纵向刚度估计可以描述为以下方程：

x(t)＝θ^T(t)k_s(t)+e(t) (26)

其中k_s(t)是随时间变化的轮胎纵向刚度；x(t)是共振频率f₀的平方；e(t)是误差函数，其取值参见(27)式；θ(t)是一个矩阵，θ^T(t)是θ(t)这个矩阵的转置，矩阵是回归系数函数，在此方法中θ(t)为是通过(20)式得到，式中R为车轮滚动半径，I为车轮转动惯量，r_x为轮胎松弛长度，在轮胎垂向力和滑移率变化较小时，可认为是常量。

对于式(26)的求解可按如下步骤进行：

步骤一：获取系统输出的x(t)，即上一模块估计共振频率的平方，并确定回归系数函数θ^T(t)。

步骤二：根据t时刻的系统真实输出值和上一时刻t-1预测的本时刻输出值来计算t时刻e(t)：

e(t)＝x(t)-θ^T(t)k_s(t-1) (27)

步骤三：利用下式求t时刻解增益向量K(t):

$K\left(t\right)=\frac{P(t-1)\mathrm{\θ}\left(t\right)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}---\left(28\right)$

式中P(t-1)为上一时刻t-1预测的本时刻的协方差矩阵，参见式(29)。χ是遗忘因子,本例中取值在(0.9，1)。

利用下式计算下一时刻t+1的协方差矩阵P(t):

$P\left(t\right)=\frac{1}{\mathrm{\χ}}[P(t-1)-\frac{P(t-1)\mathrm{\θ}\left(t\right){\mathrm{\θ}}^T\left(t\right)P(t-1)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}]---\left(29\right)$

步骤四：计算轮胎纵向刚度：

k_s(t)＝k_s(t-1)+K(t)e(t) (30)

因为轮胎纵向刚度与附着系数是线性相关的，本发明通过估计电动轮系统的共振频率，进而根据轮胎纵向刚度与电动轮共振频率之间关系式计算得到纵向刚度，进一步根据纵向刚度估计出来附着系数。

上述实施例仅用于说明本发明，其中方法实施的步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (3)

1.一种基于电机与车轮耦合特性的路面附着系数估计方法，其特征在于，包括以下步骤：

1)首先建立轮胎纵向刚度与车轮共振频率之间的关系，为：

$f_0\≈\frac{R}{2\mathrm{\π}}\sqrt{\frac{k_s}{{\mathrm{Ir}}_x}}---\left(20\right)$

式中f₀为共振频率，R为车轮滚动半径，k_s为轮胎纵向刚度，I为车轮转动惯量，r_x为轮胎的纵向松弛长度；

2)然后在整车控制器获取实时的车轮转速ω和电机的电流i_q基础上，利用MATLAB中的nonlinear ARX mdel模块，输出二阶系统模型传递函数的系数a₁，a₂，a₃，然后利用公式(22)找到两个解λ_i(1，2）：G(λ)＝a₁λ²+a₂λ+a₃ (22)

式中λ是拉布拉斯算子，a₁，a₂，a₃为系数；

3)再按照公式(23)、(24)、(25)计算共振频率f₀：

$f_0=\frac{\sqrt{d_i^2-c_i^2}}{2\mathrm{\π}}---\left(23\right)$

式中

$d_i=-\frac{1}{2\mathrm{\Δ}T}\frac{\mathrm{Im}\left({\mathrm{\λ}}_i\right)}{\mathrm{Re}\left({\mathrm{\λ}}_i\right)}---\left(25\right)$

上述ΔT为采样时间，Re、Im分别表示数学计算中的实部和虚部；

4)将公式(23)取得的f₀代入到公式(20)中，在公式(20)的基础上估计轮胎纵向刚度k_s；

5)在任意时刻，将得到的轮胎纵向刚度输入到路面附着系数估计模块中，计算得到路面附着系数：

${\mathrm{\μ}}_{max}\left(t\right)=a_{{\mathrm{\μ}}_{max}}{k}_s\left(t\right)+b_{{\mathrm{\μ}}_{max}}$

式中的系数根据实验数据确定，采用数据拟合的方法得到。

2.根据权利要求1所述的路面附着系数估计方法，其特征在于，在建立轮胎纵向刚度与车轮共振频率之间的关系时，是按照如下方法实现的：

1)建立单轮动力学模型：

$I\stackrel{\·}{\mathrm{\ω}}=T_d-{F_x}^{D}R---\left(1\right)$

式中I为车轮转动惯量，ω为车轮转速，是ω关于时间的导数，T_d为电机输出转矩，F_x ^D是车轮瞬态纵向力，R为车轮滚动半径；

2)建立稳态轮胎模型：

$F_x^S=k_{s}s+F_{x0}$

式中是车轮稳态纵向力，k_s为轮胎纵向刚度，s为滑移率，不同情况取不同值，F_x0为滑移率为零时的纵向力，v为车辆纵向速度，ω为车轮转速，R为车轮滚动半径；

3)建立瞬态轮胎模型：

$\mathrm{\τ}{{\stackrel{\·}{F}}_x}^D+{F_x}^D={F_x}^S$

$\mathrm{\τ}=\frac{r_x}{v}---\left(3\right)$

式中F_x ^D为车轮瞬态纵向力，是F_x ^D关于时间的导数，τ为时间常数，为车轮稳态纵向力，r_x为轮胎的纵向松弛长度；

4)建立简化电机模型：

T_d＝Ki_q (4)

式中i_q为电机的电流，K为比例常数，通过实验测得或电机厂商提供，T_d为电机输出转矩；

5)在上述几种车辆模型和简化电机模型的基础上，假设电机转矩由两部分组成，分为恒定部分与高频部分，如下所示：

T_d＝T₀+T₁sin(2πf·t) (5)

式中T₀表示恒定转矩，认为是一个相对恒定的值；而T₁sin(2πf·t)表示高频转矩，其中T₁为转矩振幅，f为高频转矩的频率，t是表示某一时刻；

6)结合公式(1)，将公式(5)带入到公式(1)中，单轮动力学模型可表示为：

$I\stackrel{\·}{\mathrm{\ω}}=T_0+T_{1}sin(2\mathrm{\π}f\·t)-F_x^{D}R---\left(6\right)$

对公式(6)两端同时对时间求导，得到公式(7)：

$I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\π}f\·T_{1}cos(2\mathrm{\π}f\·t)-{\stackrel{\·}{F}}_x^{D}R---\left(7\right)$

再将公式(7)乘以时间常数τ，然后与公式(6)求和得公式(8)：

$\begin{array}{c}I\mathrm{\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+T_1\sin (2\mathrm{\π}f\·t)-F_x^{D}R\\ +2\mathrm{\π}f\mathrm{\τ}\·T_1\cos (2\mathrm{\π}f\·t)-\mathrm{\τ}{\stackrel{\·}{F}}_x^{D}R\\ =T_0+T_1\sin (2\mathrm{\π}f\·t)-F_x^{S}R\\ +2\mathrm{\π}f\·T_1\cos (2\mathrm{\π}f\·t)\end{array}---\left(8\right)$

结合公式(2)、(3)和三角函数公式，公式(8)可简化为公式(9)：

$\begin{array}{c}I\mathrm{\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+{\mathrm{aT}}_1\sin (2\mathrm{\π}f\·t+\mathrm{\φ})-(k_{s}s+F_{x0})R\\ a=\sqrt{1+{\left(2\mathrm{\π}f\right)}^2},\mathrm{\φ}=\arctan \left(2\mathrm{\π}f\right)\end{array}---\left(9\right)$

式中aT₁为合并后高频信号的振幅，φ为合并后高频信号的初始相位；

对于驱动工况，考虑其滑移率定义，公式(9)可以表示为公式(10)：

$I\mathrm{\τ}\stackrel{\·\·}{\mathrm{\ω}}+I\stackrel{\·}{\mathrm{\ω}}=T_0+{\mathrm{aT}}_{1}sin(2\mathrm{\π}f\·t+\mathrm{\φ})-k_{s}R\frac{R\mathrm{\ω}-v}{R\mathrm{\ω}}-F_{x0}R---\left(10\right)$

进一步对公式(10)两端求导可得公式(11)：

$I\mathrm{\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=a2\mathrm{\π}f\·T_{1}cos(2\mathrm{\π}f\·t+\mathrm{\φ})+k_{s}R\frac{\stackrel{\·}{v}\·R\mathrm{\ω}-R\stackrel{\·}{\mathrm{\ω}}\·v}{{\left(R\mathrm{\ω}\right)}^2}---\left(11\right)$

假设Rω≈v，由于车辆的惯性远大于车轮的惯性，所以车辆的加速度相比车轮的角加速度可以忽略，公式(11)可进一步简化为公式(12)：

$I\mathrm{\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}=2\mathrm{\π}fa\·T_{1}cos(2\mathrm{\π}f\·t+\mathrm{\φ})-\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}---\left(12\right)$

令：

T₂＝T₁cos(2πf·t+φ) (13)

公式(12)可表示为公式(14)：

$I\mathrm{\τ}\stackrel{\·\·\·}{\mathrm{\ω}}+I\stackrel{\·\·}{\mathrm{\ω}}+\frac{k_s{R}^2}{v}\stackrel{\·}{\mathrm{\ω}}=2\mathrm{\π}fa\·T_2---\left(14\right)$

对公式(14)两端作拉普拉斯变换，得到公式(15)：

${\mathrm{I\τ\λ}}^3\mathrm{\ω}\left(\mathrm{\λ}\right)+{\mathrm{I\λ}}^2\mathrm{\ω}\left(\mathrm{\λ}\right)+\frac{k_s{R}^2}{v}\mathrm{\λ}\mathrm{\ω}\left(\mathrm{\λ}\right)=2{\mathrm{\πfaT}}_2\left(\mathrm{\λ}\right)---\left(15\right)$

式中λ为拉普拉斯算子；

这样就得到了电机转矩到轮速的传递函数，公式(16)：

$\frac{T_2\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{{\mathrm{I\τ\λ}}^3+{\mathrm{I\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\π}f\sqrt{1+{\left(2\mathrm{\π}f\right)}^2}}---\left(16\right)$

结合公式(4)和(16)，得到电机电流到轮速的传递函数，公式(17)：

$\frac{i_q\left(\mathrm{\λ}\right)}{\mathrm{\ω}\left(\mathrm{\λ}\right)}=\frac{{\mathrm{I\τ\λ}}^3+{\mathrm{I\λ}}^2+\frac{k_s{R}^2}{v}\mathrm{\λ}}{2\mathrm{\π}fK\sqrt{1+{\left(2\mathrm{\π}f\right)}^2}}---\left(17\right)$

令λ＝j2πf，其中j表示虚部，合并同类项，求模即可得到电机电流到轮速的幅频函数，公式(18)：

$\begin{array}{c}\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}=\left|\frac{i_q\left(j2\mathrm{\π}f\right)}{\mathrm{\ω}\left(j2\mathrm{\π}f\right)}\right|\\ =\frac{|-I{\left(2\mathrm{\π}f\right)}^2+j(\frac{k_s{R}^2}{v}2\mathrm{\π}f-I\mathrm{\τ}{\left(2\mathrm{\π}f\right)}^3)|}{2\mathrm{\π}fK\sqrt{1+{\left(2\mathrm{\π}f\right)}^2}}\\ =\frac{|-I\left(2\mathrm{\π}f\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\π}f\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\π}f\right)}^2}}\end{array}---\left(18\right)$

作如下近似：当电机车轮系统发生共振时，即有公式(19)：

$\begin{array}{c}\min \[\frac{A\left(i_q\right)}{A\left(\mathrm{\ω}\right)}\]=\min \[\frac{|-I\left(2\mathrm{\π}f\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\π}f\right)}^2}{v}\right)|}{K\sqrt{1+{\left(2\mathrm{\π}f\right)}^2}}\]\\ \≈\min \[\frac{|-I\left(2\mathrm{\π}f\right)+j\left(\frac{k_s{R}^2-{\mathrm{Ir}}_x{\left(2\mathrm{\π}f\right)}^2}{v}\right)|}{K\·2\mathrm{\π}f}\]\\ =\min \[|\frac{-I}{K}+j\left(\frac{k_s{R}^2/2\mathrm{\π}f-{\mathrm{Ir}}_x\left(2\mathrm{\π}f\right)}{vK}\right)|\]\end{array}---\left(19\right)$

应用求极值的方法，得到其最小值对应的频率，即共振频率，见公式(20)。

3.根据权利要求1或2所述的路面附着系数估计方法，其特征在于，基于递归最小二乘法进行随时间变化的轮胎纵向刚度估计，递归最小二乘法的轮胎纵向刚度估计描述为以下方程：

其中k_s(t)是随时间变化的轮胎纵向刚度；x(t)是共振频率f₀的平方；e(t)是误差函数，是矩阵的转置，矩阵是回归系数函数，为式中R为车轮滚动半径，I为车轮转动惯量，r_x为轮胎的纵向松弛长度；

对式(26)求解k_s(t)按如下步骤进行：

步骤一：获取系统输出的x(t)，并确定回归系数函数

步骤二：根据t时刻的系统真实输出值和上一时刻t-1预测的本时刻输出值来计算t时刻e(t)：

步骤三：利用下式求t时刻解增益向量K(t):

$K\left(t\right)=\frac{P(t-1)S\left(t\right)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}---\left(28\right)$

式中P(t-1)为上一时刻t-1预测的本时刻的协方差矩阵，参见式(29)，χ是遗忘因子,本例中取值在(0.9，1)，

利用下式计算下一时刻t+1的协方差矩阵P(t):

$P\left(t\right)=\frac{1}{\mathrm{\χ}}\[P(t-1)-\frac{P(t-1)\mathrm{\θ}\left(t\right){\mathrm{\θ}}^T\left(t\right)P(t-1)}{\mathrm{\χ}+\mathrm{\θ}{\left(t\right)}^{T}P(t-1)\mathrm{\θ}\left(t\right)}\]---\left(29\right)$

步骤四：计算轮胎纵向刚度：

k_s(t)＝k_s(t-1)+K(t)e(t) (30)。