顾亮, 王振宇, 王振峰
北京理工大学 机械与车辆学院, 北京 100081
收稿日期：2017-08-21
基金项目：国家自然科学基金资助项目(U1564210)；中国博士后科学基金资助项目(2016M600934，BX201600017)；国家留学基金委项目(CSC201706030029)。
作者简介：顾亮(1957-)，男，山东淄博人，北京理工大学教授，博士生导师。

摘要：针对车辆悬架状态无法准确估计的问题, 设计了自适应交互式多模型卡尔曼滤波(IMMKF)状态观测器.首先, 建立了标准路面激励模型与四分之一线性化悬架模型; 然后, 利用递归最小二乘方法与IMMKF理论, 设计了不同工况下广义悬架模型自适应IMMKF状态观测器; 最后，分析了在标准C级路面激励工况下簧载质量变化对悬架系统状态估计的影响.仿真与台架试验结果表明, 在簧载质量变化工况下, 所设计的自适应IMMKF状态观测器与传统卡尔曼滤波状态观测器相比其估计精度至少可以提高20%.
关键词：状态估计交互式多模型卡尔曼滤波递归最小二乘算法悬架系统簧载质量
Vehicle Suspension System State Estimation Combining with Interacting Multiple Model Kalman Filter
GU Liang, WANG Zhen-yu, WANG Zhen-feng
School of Mechanical and Vehicle Engineering, Beijing Institute of Technology, Beijing 100081, China
Corresponding author: WANG Zhen-feng, E-mail: wangzhenfeng44827@163.com
Abstract: In order to estimate accurately the vehicle suspension state, an interactive multiple model adaptive Kalman filter(IMMKF)state observer was proposed. Firstly, a standard road excitation model and a quarter vehicle linear model were established. Secondly, by combining recursive least square algorithm with IMMKF theory, an IMMKF state observer was designed based on the augmented suspension model in various working conditions. Finally, the influence on the state estimation of the suspension system with the change of sprung mass under the ISO level C road input excitation was analyzed. The results of simulation and experiment on a quarter of vehicle test rig showed that compared with the tradition Kalman filter(KF)state observer, the estimation accuracy of the proposed IMMKF state observer could be improved beyond 20% with the change of sprung mass.
Key words: state estimationinteracting multiple model Kalman filter(IMMKF)recursive least square algorithmsuspension systemsprung mass
利用可控悬架对车辆行驶性能进行优化控制时, 悬架系统状态的精确程度直接影响可控悬架的控制效果, 但是由于技术与成本的限制, 有些控制参数无法准确测量或传感器价格昂贵, 并且可靠性得不到保证, 难以获取所需的状态参数.基于上述情况, 需要搭建合适的观测器对悬架系统进行状态估计.考虑到典型路面干扰产生的随机激励所引起的垂向车辆动力学的随机性, 多采用卡尔曼滤波来进行车辆的状态估计^[1-2].但是传统的卡尔曼滤波器有严格的条件假设:模型应该是精确的, 统计噪声应该是高斯白噪声.因此, 仅采用在不同工作条件下的传统卡尔曼滤波器难以得到更高精度的状态估计.
目前普遍使用的多模型估计方法是基于多模型的卡尔曼滤波算法^[3]，例如:伪贝叶斯算法^[4]、交互式多模型(IMM)算法^[5]和交互式多模型自适应估计(MMAE)^[6].交互式多模型算法已经被证实可以以较低计算成本实现良好的计算效果^[7].交互式多模型算法通常将模式转移建模为与连续状态无关的恒定模式概率转移的马尔科夫过程; 在很多情况下, 当连续状态位于连续状态空间的子集内时就会发生模式转移^[8].
在基于标称系统模型进行状态观测估计时, 不确定系统参数可能造成估计质量的大幅降低.状态估计的误差就会变大.因此, 本文使用递归最小二乘法(RLS)滤波来进行参数识别, 并且这种方法已被证实是一种有效的参数识别算法^[9].文献[10]分析了一种用于参数识别的时变系统的并行自动谐调递归最小二乘算法; 文献[11]研究了一种引入遗忘因子的不稳定递归最小二乘算法与卡尔曼滤波器相结合来估计低频机电模型的观测器.
根据上述分析, 本文通过考虑不同模型之间的模式转移, 提出了一种新的基于交互式多模型卡尔曼滤波与递归最小二乘算法耦合的状态估计方法.
1 道路激励模型本研究假定路面不平度是一种均匀的各向同性的高斯随机过程, 它的统计特征可以用功率谱密度(PSD)来描述.根据ISO 8601^[12], 路面功率谱密度定义为

(1)

其中:n(单位为m^-1)为空间频率, n₀=0.1 m^-1是参考空间频率; G_q(n₀)(单位为m^-3)为参考频率n₀下的路面功率谱密度值, 称为路面不平度系数; w为频率指数, 它决定路面功率谱密度的频率结构, 通常取w=2^[13].
考虑到角频率ω, 式(1)可写成

(2)

基于频谱分解定理, G(ω)在功率谱密度函数稳定随机控制过程中起作用.式(2)因此可以表示成式(3)和式(4).

(3)

(4)

根据路面生成的时间频率公式^[14]：

(5)

其中σ是白噪声信号的系数.
式(4)和式(5)可以表示成式(6)和式(7).

(6)

(7)

, 由式(6)得出：

(8)

2 车辆悬架模型2.1 四分之一车辆模型四分之一线性车辆模型如图 1所示.
图 1(Fig. 1)

图 1 四分之一车辆悬架模型Fig.1 Model of quarter vehicle suspension

模型动力学方程可表示为

(9)

其中:F_d表示可控阻尼力; k_s, k_t表示悬架弹簧刚度和轮胎弹簧刚度; x_b, 和分别是簧上质量的位移、速度和加速度; x_w, 和分别是簧下质量的位移、速度和加速度; x_r是道路不平度, 悬架系统的振动由路面激励产生.
系统的状态向量和输出向量为
其中, 4种状态分别是频响空间、轮胎侧偏、簧上质量速度和簧下质量速度.
因此, 状态空间变量的状态空间方程可写为式(10)：

(10)

其中:
其中:c表示悬架阻尼; w和v表示假定独立的过程高斯噪声和观测高斯噪声; 代表路面不平度的速度变化.Q(k)=E(ww^T)和R(k)=E(vv^T)分别是相应过程噪声方差和观测噪声协方差.
状态空间离散表达式为

(11)

式中:k表示离散时间的采样时刻; T表示采样时间.
2.2 广义车辆悬架模型分析悬架的非线性特性主要由弹簧特性、递减非对称的悬架阻尼力：速度的特性、库仑摩擦和悬架支点与车轮平面的运动传递系数产生.这些非线性特性可由相关的非线性阻尼c(k)和变化的弹簧质量mb(k)导出.基于图 1采用一系列不同的附加载荷可以推导出车辆悬架的簧上质量运动方程m_b(k)=m_b(k-1)+Δm(k).
本研究使用马尔科夫线性四分之一车辆模型.因此, 广义的车辆四分之一状态矢量空间模型可由式(12)表示.

(12)

其中:
其中, f_1i和f_2i分别表示第i个不变的最小区域空间和最大区域空间的组合.对于路面速度来说, 未知的干扰可以用高斯白噪声过程来近似.模型的不确定性部分也可以假设成高斯白噪声.每一个线性子模型i都可依据表 1所示的数值结合状态空间表达式联立表示.
表 1(Table 1)

表 1 四分之一悬架参数Table 1 Parameters for the quarter suspension 悬架参数 符号/单位 数值 簧载质量 mb/kg 420 非簧载质量 mw/kg 39 悬架弹簧刚度 Ks/(N·m-1) 20 000 轮胎刚度 Kt/(N·m-1) 185 000 最大阻尼 c/(N·s·m-1) 2 000 总簧载质量 m/kg 450~650

表 1 四分之一悬架参数 Table 1 Parameters for the quarter suspension

3 递归最小二乘算法和交互式多模型卡尔曼滤波算法分析3.1 簧载质量的递归最小二乘算法递归最小二乘算法是根据已知数据来估计系统参数.假设系统为

(13)

其中:Θ表示输入矩阵; Y表示输出矩阵; θ表示待识别参数矩阵.
本文提出了可以变换为递归最小二乘法的应用程序.k时刻的簧上质量可以写为

(14)

其中, Θ(k)定义为在第k步中可测量的数据因子.y(k), m_b(k)和e(k)分别表示第k步中输出的测量值、被识别的参数和高斯白噪声过程.
基于上述分析, k时刻簧上质量的递归格式可写为

(15)

联立更新的观测变量和误差协方差矩阵R, 最小二乘公式可写为

(16)

(17)

其中, K_k(k)代表增益矩阵.
修正了它的估计.图 2给出递归最小二乘算法流程图.
图 2(Fig. 2)

图 2 递归最小二乘算法流程图Fig.2 Flow chart of RLS algorithm

3.2 交互式多模型卡尔曼滤波算法理论的表示本文提出了一个时间离散随机线性混合系统, 如式(18)所示:

(18)

其中:k=1, 2, …; x(k)∈Rⁿ和z(k)∈R^p分别是离散状态向量和观测向量; m(k)∈{1, 2, …, r}是离散状态(或模式); A_m(k)和C_m(k)是与k时刻的m(k)一致的状态转移矩阵^[13].过程噪声w_m(k)(k)和观测噪声v_m(k)(k)是互不相关的零均值高斯白噪声.高斯序列:Q(k)= E(w_m(k)w_m(k)^T)和R(k)= E(v_m(k)v_m(k)^T)分别是它们相应的过程噪声方差矩阵和观测噪声协方差矩阵.模式转移即m(1), m(2), …, m(k)是一个马尔科夫链, 可以写成一个连续态非独立的模态转移矩阵:

(19)

其中，π_ij(x(k))是一个在连续态x(k)的条件下从i模态到j模态的条件概率.

(20)

本研究用P[·|·]来表示条件概率密度分布函数(PDF).
所以, 交互式多模型卡尔曼滤波方法步骤如下:
第一步, 为j模态相同等级的模态定义交互式输入.
由前一步的目标状态估计和概率模型(μ_j(k-1))求得当前步骤下悬架系统的混合模态和协方差估计, P_oj(k-1|k-1).
预测模态j的概率:

(21)

从i模态到j模态的混合概率:

(22)

j模态的混合状态估计:

(23)

j模态的混合协方差估计:

(24)

其中：P_ij是从i模态到j模态的转移概率；μ_j(k-1)是j模态在t=k-1步的概率.
第二步, 计算j模态的卡尔曼滤波.
为了更新悬架的预测状态估计(, P_j(k|k)), 本研究使用卡尔曼滤波来输入状态, P_oj(k-1|k-1)和z(k).
预测从k-1时刻到k时刻的状态改变:

(25)

其中A_j为系统离散矩阵.
预测从k-1时刻到k时刻的协方差改变:

(26)

计算k时刻卡尔曼滤波的增益方程:

(27)

计算卡尔曼滤波方程:

(28)

计算卡尔曼滤波器的协方差:

(29)

第三步, 更新j模态的概率.
采用似然函数来更新模型的概率μ_j(k).j模态的表达式可改写为

(30)

在此基础上, j模态的概率:

(31)

第四步, 定义j模态相同等级的交互式输出.
基于j模态的概率, 整体的模态估计可通过加权每一个滤波器的状态估计结果得到.
计算整体的模态估计:

(32)

计算整体的协方差估计:

(33)

交互式多模型卡尔曼滤波器的流程图如图 3所示.
图 3(Fig. 3)

图 3 交互式多模型卡尔曼滤波算法的流程图Fig.3 Flow chart of IMMKF algorithm

4 仿真和试验4.1 仿真结果在本次仿真中, 用第3节推导出的交互式多模型卡尔曼滤波算法和递归最小二乘算法耦合的状态估计方案来估计悬架系统垂向动力学的广义状态.本研究假设条件为:轮胎不会与地面脱离接触; 并且假定马尔科夫转移矩阵在仿真之前已知; 道路激励是30 km/h的ISO-C级谱.
本研究在仿真中比较上文提到的状态参数交互式多模型卡尔曼滤波联合算法和单一状态的卡尔曼滤波器状态估计性能.计算出状态估计误差的均方根值, 并在表 2中总结了两种状态估计方法的仿真结果.通过比较仿真的状态估计精度, 所有相关信号的性能的变化即均方根误差波动始终保持在19%内.
表 2(Table 2)

表 2 以30 km / h的速度在多种道路C级不平度的条件下不同的卡尔曼滤波均方根值误差估计Table 2 Estimation RMS error of different KF various on road level C profile at velocity of 30 km/h 滤波器状态估计 滤波方式 参考仿真结果均方根误差/% xb-xw KF 7.78 IMMKF 0.48 xw-xr KF 18.14 IMMKF 3.95 KF 6.48 IMMKF 0.65 KF 0.64 IMMKF 0.31

表 2 以30 km / h的速度在多种道路C级不平度的条件下不同的卡尔曼滤波均方根值误差估计 Table 2 Estimation RMS error of different KF various on road level C profile at velocity of 30 km/h

研究了增加Δm_b=200 kg的额外车体负载的情况.这一载荷近似增加了约45%的车体质量, 可以看作是在汽车上坐了8名乘客.图 4阐明了卡尔曼滤波和交互式多模型卡尔曼滤波的相应结果.如表 2和图 4所示, 交互式多模型卡尔曼滤波可以提高多模型悬架的状态估计精度.
图 4(Fig. 4)

图 4 道路等级C(速度30 km / h)的结果Fig.4 Results of road level C(at velocity of 30 km/h) (a)—xb-xw; (b)—xw-xr; (c)—; (d)—; (e)—error(xb-xw); (f)—error(xw-xr); (g)—; (h)—.

4.2 试验验证卡尔曼滤波和交互式多模型卡尔曼滤波算法性能估计在图 5所示的四分之一车辆悬架的试验台架上进行.在状态估计的试验过程中, 道路激励由一个液压油缸施加在车轮上, 簧上质量和簧下质量的加速度数据来自道路等级ISO-C级激励谱.此试验中, 数量估计的道路激励参考信号是使用在试验设备中的第二部分通过离线计算得到的.
图 5(Fig. 5)

图 5 四分之一车辆悬架试验台架状态估计Fig.5 Quarter vehicle suspension test rig for state estimation

4.2.1 道路模拟机参数道路模拟机的参数如表 3所示.
表 3(Table 3)

表 3 道路模拟机参数Table 3 Parameters of road simulation machine 仪器性能指标 仪器参数 试验台架最大激励质量/kg 1 000 公称激励力/kN 160 公称激励位移/mm ±200 最大速度/(m·s-1) 2.5 最大加速度/g 16 工作频率/Hz 0.1~50

表 3 道路模拟机参数 Table 3 Parameters of road simulation machine

4.2.2 试验结果分析在状态估计中, 仿真结果均方根误差和试验数据均方根误差如表 4所示.在表 4中, 仿真结果与试验结果吻合良好, 可见仿真和试验结果之间的最大均方根误差小于21%, 位于合理范围内.
表 4(Table 4)

表 4 试验数据和仿真数据Table 4 Test data and simulation data 滤波器状态估计 滤波方式 xb-xw KF 7.78 8.59 IMMKF 0.48 0.65 xw-xr KF 18.14 20.80 IMMKF 3.95 5.89 KF 6.48 8.96 IMMKF 0.65 0.97 KF 0.64 0.83 IMMKF 0.31 0.37

表 4 试验数据和仿真数据 Table 4 Test data and simulation data

比较试验结果和仿真结果可见, 使用本文提出的算法可以提高悬架系统的状态估计精度.
5 结论1) 利用交互式多模型卡尔曼滤波算法与递归最小二乘算法耦合得到的自适应卡尔曼滤波方法, 可以在有效识别悬架参数的基础上, 精确得到悬架系统状态.仿真结果说明, 通过与传统卡尔曼滤波方法的比较, 所提出的新方法能够获得更高精度的悬架系统状态估计.
2) 仿真与试验结果表明, 在车辆四分之一悬架台架试验中, 当悬架簧载质量发生变化时, 传统卡尔曼滤波方法估计性能会大幅恶化, 所以此工况下可以利用本文所提出的新方法进行悬架实际状态的实时精确估计.通过仿真与试验可知, 此方法的仿真和试验结果之间的状态估计均方根误差小于21%.
参考文献

[1]	Yu F, Crolla D. An optimal self-tuning controller for an active suspension[J].Vehicle System Dynamics, 1998, 29: 51–65.DOI:10.1080/00423119808969366
[2]	Rigatos G, Siano P, Pessolano S. Design of active suspension control system with the use of Kalman filter-based disturbances estimator[J].Cybernetics and Physics, 2012, 1: 259–275.
[3]	Li X R, Jilkov V. Survey of maneuvering target tracking.part Ⅴ:multiple-model methods[J].IEEE Transactions on Aerospace and Electronic Systems, 2005, 41: 1255–1321.DOI:10.1109/TAES.2005.1561886
[4]	Maybeck P. Stochastic models, estimation, and control[M]. New York: Academic Press, 1982.
[5]	Blom H, Bar-Shalom Y. The interacting multiple model algorithm for systems with Markovian switching coefficients[J].IEEE Transactions on Automatic Control, 1988, 33: 780–783.DOI:10.1109/9.1299
[6]	Hanlon P, Maybeck P. Multiple-model adaptive estimation using a residual correlation Kalman filter bank[J].IEEE Transactions on Aerospace and Electronic Systems, 2000, 36: 393–406.DOI:10.1109/7.845216
[7]	Bar-Shalom Y, Kirubarajan X R, Li T. Estimation with applications to tracking and navigation[M]. New York: Wiley, 2001.
[8]	Campo L, Mookerjee P, Bar-Shalom Y. State estimation for systems with sojourn-time-dependent Markov model switching[J].IEEE Transactions on Automatic Control, 1991, 36: 238–243.DOI:10.1109/9.67304
[9]	Zhu Y. Efficient recursive state estimator for dynamic systems without knowledge of noise covariance[J].Mechanic Engineering, 2003, 8: 64–66.
[10]	Li J S, Zheng Y J, Lin Z P. Recursive identification of time-varying systems:self-tuning and matrix RLS algorithms[J].Systems & Control Letters, 2014, 66: 104–110.
[11]	Moreno I, Messina A R. Adaptive tracking of system oscillatory modes using an extended RLS algorithm[J].Electric Power Systems Research, 2014, 114: 28–38.DOI:10.1016/j.epsr.2014.03.034
[12]	ISO.Mechanical vibration-road surface profiles-reporting of measured data: ISO8601-1995[S].Melbourne: ARRB Group Ltd., 1995.
[13]	余志生. 汽车理论[M]. 北京: 机械工业出版社, 2009: 207-212. ( Yu Zhi-sheng. Auto theory[M]. Beijing: China Machine Press, 2009: 207-212.)
[14]	Zhang G, Fang Z, Chen S, et al. White noise simulation for road roughness based on power function[J].Automotive Engineering, 2008, 30: 44–47.