履带车辆被广泛地应用于行星探索、采矿、军事以及农业等领域^[1-5]。由于采用差速转向机构，履带车辆具有优越的机动性和通过性，对复杂的非结构化道路具有很强的地形适应能力。但是，履带接地面积较大，由于滑动转向产生的相互作用力很复杂，为履带车辆建立精确的模型较为困难。对于智能履带车辆，无论是路径规划还是横向控制都十分依赖精确的履带车辆模型^[6-7]。因此，进一步提升履带车辆模型精度具有重要意义和应用价值。

目前，许多学者基于滑动参数估计的方法来进一步研究履带车辆模型^[8-9]。现有的方法可分为2类。第1类研究集中于基于履带与地面间相互作用力的动力学模型，其核心思想是将履带与地面间的作用力表示为滑动参数的函数。Ward等^[10]针对差动转向的轮式车辆，将轮胎牵引力与制动力引入动力学模型，并结合惯性测量单元与全球卫星定位系统估计车辆的滑动参数。Ray等^[11]提出了一种将半经验的地形力学模型与车辆履带的相互作用进行净效应解耦，并通过采集传感器数据和离线处理数据来识别履带滑移特性的方法。Dar等^[12]利用关键土壤参数推导的牵引力模型，通过带有状态噪声补偿的扩展Kalman滤波器(extended Kalman filter, EKF)来估计滑动参数。周波等^[13]建立了基于运动学与动力学的修正模型，利用无迹Kalman滤波(unscented Kalman filter, UKF)对履带车辆的位姿与滑动参数进行联合估计。但是，由于地面与履带复杂的力学关系具有高度的非线性，因此基于动力学模型的参数估计方法几乎无法应用于对实时性能需求很高的智能履带车辆系统。

第2类研究主要关注基于履带车辆运动学模型的滑动参数估计。Song等^[14]利用滑模观测器与EKF观测器相结合的计算机视觉对车辆滑移率进行估计。为了提高基于模型预测控制的路径跟踪精度，Burke采用自适应最小二乘估计了车辆滑移率^[15]。焦俊等^[16]针对农用履带车辆，提出了基于UKF的滑移率在线估计方法。但是，若采用基于滑移率的运动学模型，当低速侧履带主动轮转速为0时，传统的运动学模型不能根据主动轮转速反求车体移动速度，这为实际应用带来了很大的局限性^[3]。

为了克服基于滑移率的运动学模型带来的局限性，Martinez等^[17]使用了基于瞬时转向中心(instantaneous center of rotation, ICR)的运动学模型，通过对履带接地段的运动分析，得到车体的侧、纵向速度及角速度不仅与两侧履带卷绕速度有关，还与履带及车体的ICR位置密切相关的结论，但是其ICR位置是使用遗传算法离线学习得到的。Pentzer等^[18]利用EKF实现了对ICR位置的实时估计，并分别利用轮式车辆与履带车辆验证了其算法的有效性。熊光明等^[19]假定ICR位置可以表示为侧向加速度和曲率绝对值的线性组合，并基于车辆的预测姿态变化和实际姿态变化之间的残差，分别应用EKF与Levenberg-Marquardt(LM)方法在线学习滑动参数，最后通过估计的ICR位置与相应的运动学模型进行轨迹预测。实验分析表明，LM方法的参数估计精度更高，而EKF的计算速度较快。赵梓烨等^[20]考虑了车辆的横向、纵向以及垂向运动，采用EKF对ICR位置进行估计，利用运动学模型构建了模型预测控制横向控制器，进一步提升了车辆的控制精度。

由于履带车辆工作条件复杂、地形随机性强等不确定性因素，履带车辆系统具有较强的非线性^[21]。传统的Kalman滤波器在很大程度上依赖于精确的车辆模型并且假设过程噪声和测量噪声为Gauss白噪声。在实际应用中，简化模型和噪声协方差不确定性可能存在时变和非线性误差，而这些误差实际上不能被认为是Gauss白噪声，并且会造成较大的估计误差。为克服简化模型和噪声协方差不确定性带来的问题，本文提出了一种基于轨迹预测补偿的履带车辆双重无迹Kalman滤波(dual-layer unscented Kalman filter, DUKF)滑动参数实时估计算法，将ICR位置作为滑动参数，并提出利用轨迹预测残差对初步估计的滑动参数进行补偿。本文比较了不同预测时域对DUKF估计精度的影响，并将结果与传统的EKF与UKF估计方法进行了对比。仿真结果表明，DUKF进一步提升了滑动参数估计的精度。除此之外，DUKF可以根据履带车辆的当前状态预测滑动参数的变化趋势，进一步提高算法的收敛速度，提升了估计算法在极限工况下的鲁棒性。

1 基于ICR的履带车辆建模

履带车辆运动过程中，履带与地面之间会产生大量的滑转滑移运动，这使得包含大量土壤参数的履带车辆动力学模型十分复杂，不适用于对实时性要求较高的路径跟踪控制与轨迹规划研究。本文采用基于ICR的履带车辆模型，如图 1所示，通过估计ICR位置就可以准确描述履带车辆滑转滑移运动^[18]。

图 1中：XOY为大地坐标系，车辆的位置与姿态定义为(x, y, φ)。xoy为车体坐标系，假设车辆的重心位于车身的几何中心。(v_x, v_y, ω)分别为车辆的纵向速度、横向速度与横摆角速度，(x_c, y_c)为车体瞬时旋转中心，(x_r, y_r)和(x_l, y_l)分别为高速侧和低速侧履带接地段的ICR，A和B为履带接地点，v_r和v_l分别为右、左侧履带接地段的纵向速度，v_A和v_B分别为A点和B点处牵连速度，L为两侧履带接地点间的距离。

在转向时，履带车辆的一些重要运动学参数大小与两侧履带和车体的ICR有关。履带接地段在任一时刻的平面运动，可看作绕ICR的旋转运动，因此A、B点的绝对速度V_A、V_B可表示为

$ \left\{\begin{array}{l} V_{A}=\left(y_{\mathrm{l}}-\frac{L}{2}\right) \omega, \\ V_{B}=\left(y_{\mathrm{r}}-\left(-\frac{L}{2}\right)\right) \omega . \end{array}\right. $

(1)

A、B点的绝对速度是牵连速度和相对速度的矢量和，A、B点的绝对速度的标量形式为

$ \left\{\begin{array}{l} V_{A}=v_{A}-v_{\mathrm{l}}, \\ V_{B}=v_{B}-v_{\mathrm{r}} . \end{array}\right. $

(2)

由于车辆作绕ICR(x_c, y_c)的旋转运动，车体上A、B的牵连速度为

$ \left\{\begin{array}{l} v_{A}=v_{x}-\frac{L}{2} v, \\ v_{B}=v_{x}+\frac{L}{2} v. \end{array}\right. $

(3)

由式(1)—(3)的履带车辆运动学关系，车身与地面间的ICR的纵坐标值为

$ y_{\mathrm{c}}=\frac{v_{x}}{\omega}. $

(4)

两侧履带与地面间的ICR纵坐标值为：

$ y_{\mathrm{l}}=-\frac{v_{\mathrm{l}}-v_{x}}{\omega}, $

(5)

$ y_{\mathrm{r}}=-\frac{v_{\mathrm{r}}-v_{x}}{\omega} . $

(6)

关于履带车辆在x轴的ICR可表示为

$ x_{\mathrm{l}}=x_{\mathrm{r}}=x_{\mathrm{c}}=-\frac{v_{y}}{\omega} . $

(7)

式(4)—(7)描述了ICR位置与履带接地点速度的运动学关系。

通过对履带接地段的运动分析，车体的横、纵向速度及横摆角速度不仅与两侧履带卷绕速度有关，还与履带及车体的ICR位置密切相关。由式(5)—(7)可以得到车体坐标系下车体纵向速度、横向速度和横摆角速度与ICR坐标值之间的关系，

$ \left\{\begin{array}{l} v_{x}=\frac{y_{\mathrm{l}} v_{\mathrm{r}}-y_{\mathrm{r}} v_{\mathrm{l}}}{y_{\mathrm{l}}-y_{\mathrm{r}}}, \\ v_{y}=\frac{v_{\mathrm{l}}-v_{\mathrm{r}}}{y_{\mathrm{l}}-y_{\mathrm{r}}} x_{\mathrm{c}}, \\ \omega=-\frac{v_{\mathrm{l}}-v_{\mathrm{r}}}{y_{\mathrm{l}}-y_{\mathrm{r}}} . \end{array}\right. $

(8)

于是，基于ICR的车辆运动微分方程为

$ \left[\begin{array}{c} \dot{x} \\ \dot{y} \\ \dot{\varphi} \end{array}\right]=\left[\begin{array}{ccc} \cos \varphi & -\sin \varphi & 0 \\ \sin \varphi & \cos \varphi & 0 \\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{c} \frac{y_{\mathrm{l}} v_{\mathrm{r}}-y_{\mathrm{r}} v_{\mathrm{l}}}{y_{\mathrm{l}}-y_{\mathrm{r}}} \\ \frac{v_{\mathrm{l}}-v_{\mathrm{r}}}{y_{\mathrm{l}}-y_{\mathrm{r}}} x_{\mathrm{c}} \\ -\frac{v_{\mathrm{l}}-v_{\mathrm{r}}}{y_{\mathrm{l}}-y_{\mathrm{r}}} \end{array}\right] . $

(9)

2 基于轨迹预测校正的DUKF滑动参数估计 2.1 无迹Kalman滤波算法

UKF是一种常用的非线性滤波算法，其核心思想是利用无迹变换(unscented transformation, UT)逼近非线性系统状态的后验概率密度函数，并通过预测和更新进行系统估计。该算法通过选取Sigma点，将这些点代入非线性方程中，重构出新的均值和方差，再将变换得到的均值、方差引入滤波器的递推过程中。

针对非线性离散系统的状态方程和观测方程，UKF递推算法如下:

$ \left\{\begin{array}{l} \boldsymbol{X}_{k+1}=f\left(\boldsymbol{X}_{k}\right)+\boldsymbol{w}_{k}, \\ \boldsymbol{Z}_{k}=h\left(\boldsymbol{X}_{k}\right)+\boldsymbol{v}_{k} . \end{array}\right. $

(10)

$ \left\{\begin{array}{l} \boldsymbol{R}=E\left(\boldsymbol{v}_{k} \boldsymbol{v}_{k}^{\mathrm{T}}\right),\\ \boldsymbol{Q}=E\left(\boldsymbol{w}_{k} \boldsymbol{w}_{k}^{\mathrm{T}}\right). \end{array}\right. $

(11)

式(10)与(11)中：X_k∈Rⁿ和Z_k∈R^m分别为状态向量与观测向量；w_k和v_k分别为被假设成Gauss白噪声的过程噪声和测量噪声; f(·)和h(·)分别为描述状态量和观测量的非线性函数，其噪声协方差矩阵分别定义为R和Q；E(·)表示括号内变量的均值。

1) 系统状态初始化。

假定履带车初始状态X₀为符合Gauss分布的随机向量，则初始状态和方差分别为

$ \left\{\begin{array}{l} \hat{\boldsymbol{X}}_{0}=E\left(\boldsymbol{X}_{0}\right), \\ \boldsymbol{P}_{0}=\operatorname{cov}\left(\boldsymbol{X}_{0}, \boldsymbol{X}_{0}^{\mathrm{T}}\right)=E\left[\left(\boldsymbol{X}_{0}-\hat{\boldsymbol{X}}_{0}\right)\left(\boldsymbol{X}_{0}-\hat{\boldsymbol{X}}_{0}\right)^{\mathrm{T}}\right] . \end{array}\right. $

(12)

其中cov(·)表示求括号内两个变量的协方差。

2) 计算加权系数。

$ \left\{\begin{array}{l} \omega_{0}^{\mathrm{m}}=\lambda /(n+\lambda), \\ \omega_{0}^{\mathrm{c}}=\lambda /(n+\lambda)+\left(1-\alpha^{2}+\beta\right), \\ \omega_{i}^{\mathrm{c}}=\omega_{i}^{m}=\frac{1}{2(n+\lambda)}, i=1,2, \cdots, 2 n . \end{array}\right. $

(13)

式(13)中：λ=α²(n+κ)-n，α和κ是确定Sigma点分布的比例参数，α的取值范围一般为1×e^-4≤α≤1，κ的取值没有特定的界限。β是一个非负的常数，对于Gauss分布来说，β=2是最优的。

3) 计算Sigma点并进行更新。

$ \left\{\begin{array}{l} \boldsymbol{\xi}_{0, k}=\hat{\boldsymbol{X}}_{k}, \quad i=0 ; \\ \boldsymbol{\xi}_{i, k}=\hat{\boldsymbol{X}}_{k}+\left(\sqrt{(n+\lambda) \boldsymbol{P}_{k}}\right), i=1,2, \cdots, n ; \\ \boldsymbol{\xi}_{i, k}=\hat{\boldsymbol{X}}_{k}-\left(\sqrt{(n+\lambda) \boldsymbol{P}_{k}}\right), i=n+1, n+2, \cdots, 2 n . \end{array}\right. $

(14)

$ \hat{\boldsymbol{X}}_{k+1 / k}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{m}} \boldsymbol{\xi}_{i, k+1 / k}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{m}} f\left(\boldsymbol{\xi}_{i, k}\right), $

(15)

$ \begin{gathered} \boldsymbol{P}_{k+1 / k}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{c}}\left(\boldsymbol{\xi}_{i, k+1 / k}-\hat{\boldsymbol{X}}_{k+1 / k}\right) \cdot \\ \left(\boldsymbol{\xi}_{i, k+1 / k}-\hat{\boldsymbol{X}}_{k+1 / k}\right)^{\mathrm{T}}+\boldsymbol{Q}_{k} . \end{gathered} $

(16)

$ \left\{\begin{array}{l} \hat{\boldsymbol{Z}}_{k+1 / k}=h\left(\hat{\boldsymbol{X}}_{k+1 / k}\right), \\ \hat{\boldsymbol{z}}_{k+1 / k}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{m}} \hat{\boldsymbol{Z}}_{k+1 / k}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{m}} h\left(\hat{\boldsymbol{X}}_{k+1 / k}\right) . \end{array}\right. $

(17)

4) 更新观测量协方差。

$ \begin{gathered} \boldsymbol{P}_{\hat{\boldsymbol{Z}}_{k+1 / k}}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{c}}\left(\hat{\boldsymbol{z}}_{k+1 / k}-\hat{\boldsymbol{Z}}_{k+1 / k}\right) \cdot \\ \left(\hat{\boldsymbol{z}}_{k+1 / k}-\hat{\boldsymbol{Z}}_{k+1 / k}\right)^{\mathrm{T}}+\boldsymbol{R}_{k} . \end{gathered} $

(18)

$ \begin{gathered} \boldsymbol{P}_{\hat{\boldsymbol{X}}_{k+1 / k} \hat{\boldsymbol{Z}}_{k+1 / k}}=\sum\limits_{i=0}^{2 n} \omega_{i}^{\mathrm{c}}\left(\boldsymbol{\xi}_{i, k+1 / k}-\hat{\boldsymbol{X}}_{k+1 / k}\right) \cdot \\ \left(\hat{\boldsymbol{z}}_{k+1 / k}-\hat{\boldsymbol{Z}}_{k+1 / k}\right)^{\mathrm{T}} . \end{gathered} $

(19)

5) 计算无迹Kalman增益并进行状态更新。

$ \boldsymbol{K}_{k+1}=\boldsymbol{P}_{\hat{\boldsymbol{X}}_{k+1 / k} \hat{\boldsymbol{Z}}_{k+1 / k}} \boldsymbol{P}_{\hat{\boldsymbol{Z}}_{k+1 / k}}^{-1}. $

(20)

$ \left\{\begin{array}{l} \hat{\boldsymbol{X}}_{k+1}=\hat{\boldsymbol{X}}_{k+1 / k}+\boldsymbol{K}_{k+1}\left(\boldsymbol{Z}_{k}-\hat{\boldsymbol{Z}}_{k+1 / k}\right), \\ \boldsymbol{P}_{k+1}=\boldsymbol{P}_{k+1 / k}-\boldsymbol{K}_{k+1} \boldsymbol{P}_{\hat{\boldsymbol{Z}}_{k+1 / k}} \boldsymbol{K}_{k+1}^{\mathrm{T}} . \end{array}\right. $

(21)

2.2 双重无迹Kalman滤波算法

由于简化模型和噪声协方差的不确定性，UKF方法在实际应用中的估计误差是不可避免的，并且在极限工况下可能会出现较大的估计误差。为了提高滑动参数的估计精度，本文提出了基于轨迹预测补偿的DUKF方法。

该方法仅利用车辆传感器信息就能在路面先验条件未知的情况下估计履带滑动参数。DUKF利用车辆的历史轨迹信息，通过上层UKF对滑动参数进行初步估计，将初步估计的滑动参数再次输入车辆模型进行轨迹预测，同时利用车辆当前的绝对速度对未来时刻的车辆轨迹进行预测，将这2种预测轨迹分别作为预测量和观测量输入下层UKF进行滑动参数的估计补偿。履带车辆滑动参数DUKF估计整体逻辑框图如图 2所示。

对于上层UKF，令滑动参数R^h=(y_l^h, y_r^h, x_c^h)^T作为非线性履带车辆系统的状态向量，其中上标h代表此变量属于上层UKF，系统输入为u^h=(v_l^h, v_r^h, φ^h)^T，观测量为Z^h=(Δx, Δy)，预测量为Ẑ^h=(Δx^h, Δy^h)。系统输入u^h为车辆历史状态信息，观测量Z^h为每个采样间隔车辆历史位姿的变化量; 预测量Ẑ^h为每个采样间隔将已选取的Sigma点、车辆历史状态信息与滑动参数R^h输入车辆运动学模型预测得到的历史相对位姿变化量，预测量Ẑ^h可通过式(22)计算：

$ \left\{\begin{array}{l} \Delta x^{\mathrm{h}}=V_{x}^{\mathrm{h}} \Delta t \cos \theta-V_{y}^{\mathrm{h}} \Delta t \sin \theta ,\\ \Delta y^{\mathrm{h}}=V_{x}^{\mathrm{h}} \Delta t \sin \theta+V_{y}^{\mathrm{h}} \Delta t \cos \theta. \end{array}\right. $

(22)

式(22)中：Δt为系统采样时间，V_x^h与V_y^h可通过式(8)计算得到。利用UKF算法迭代处理预测量Ẑ^h、观测量Z^h与非线性履带车辆系统的状态向量R^h，通过上层UKF迭代运算初步估计履带车辆滑动参数$\hat{\boldsymbol{R}}^{h}$。

对于下层UKF，令滑动参数R^p=(y_l^p, y_r^p, x_c^p)^T作为非线性履带车辆系统的状态向量，其中上标p代表此变量属于下层UKF，迭代初值R₀^p=$\hat{\boldsymbol{R}}^{\mathrm{h}}=\left(y_{1}^{\mathrm{h}}, y_{\mathrm{r}}^{\mathrm{h}}, x_{\mathrm{c}}^{\mathrm{h}}\right)$。系统输入为车辆当前时刻的状态u^p=(v_l^p, v_r^p, φ^p)^T，定义观测量为Z_k^p=(Δx_v^p, Δy_v^p)，Z_k^p利用车辆当前时刻的绝对速度与横摆角速度来预测未来一段采样间隔内车辆位姿的变化量，观测量Z_k^p的计算可表示为

$ \left\{\begin{array}{l} \Delta x_{v}^{\mathrm{p}}=V \Delta t \cos \theta_{k}, \\ \Delta y_{v}^{\mathrm{p}}=V \Delta t \sin \theta_{k}, \\ \theta_{k}=\omega_{v} k \Delta t+\theta_{0} . \end{array}\right. $

(23)

式(23)中：V、ω_v分别为传感器测得的当前时刻的车辆绝对速度与角速度，θ₀为车辆初始横摆角，θ_k为预测的第k个采样时刻车辆横摆角。

定义预测量为Ẑ_k^p=(Δx_R^p, Δy_R^p)。Ẑ_k^p为将已选取的Sigma点、车辆当前时刻的状态信息u^p与滑动参数R^p输入车辆运动学模型得到的未来一段采样间隔内相对位姿变化量，预测量Ẑ_k^p的计算可表示为

$ \left\{\begin{array}{l} \Delta x_{\boldsymbol{R}}^{\mathrm{p}}=V_{x}^{\mathrm{p}} \Delta t \cos \theta_{k}-V_{y}^{\mathrm{p}} \Delta t \sin \theta_{k}, \\ \Delta y_{\boldsymbol{R}}^{\mathrm{p}}=V_{x}^{\mathrm{p}} \Delta t \sin \theta_{k}+V_{y}^{\mathrm{p}} \Delta t \cos \theta_{k} . \end{array}\right. $

(24)

式(24)中，V_x^p与V_y^p可通过式(8)计算。利用下层UKF算法迭代处理预测量Ẑ_k^p、观测量Z_k^p与非线性履带车辆系统的状态向量R^p，通过迭代运算估计补偿后的履带车辆滑动参数ŷ_l^p、ŷ_r^p、$\hat{x}_{\mathrm{c}}^{\mathrm{p}}$，最终系统输出履带车辆滑动参数$\hat{\boldsymbol{R}}^{\mathrm{p}}$。

3 仿真验证与分析

为了验证本文所提出的DUKF滑动参数估计方法的有效性与实用性，基于MATLAB/Simulink与RecurDyn搭建了仿真模型。地形设置为平地黏土。RecurDyn模型如图 3所示，履带车辆参数如表 1所示。

由于滑动参数估计的目的是使当前状态下车辆模型更加准确，为了比较所提出方法的估计精度，本文利用结合滑动参数估计的履带车辆模型进行前向轨迹预测，预测时间为2 s，然后将预测的轨迹与2 s后车辆实际的轨迹进行对比并计算残差，即两者的残差反映了参数估计的精度。为了充分验证DUKF估计滑动参数的效果，本文除了将DUKF与EKF和UKF进行对比之外，还分别比较了将下层UKF中预测时域参数分别设置为0.1、0.5与1 s情况下的滑动参数估计结果，以说明下层UKF中不同预测时域参数对滑动参数估计精度的影响。

第1组参考轨迹由1段直线和1段圆弧组成，路径起点的坐标为(0, 0)。路径为小曲率工况，其曲率在10 s时存在突变，如图 4所示。履带车辆两侧履带卷绕速度如图 5所示，高速侧履带卷绕速度为4.8 m/s，低速侧履带卷绕速度为4.5 m/s。

DUKF的滑动参数估计结果如图 6所示，在车辆刚起步时，y_l与y_r均处于两侧履带中心，之后随着车速的增加，y_l与y_r均向履带外侧移动，滑动参数$\hat{\boldsymbol{R}}=(0.9 \mathrm{~m}, -0.9 \mathrm{~m}, 0 \mathrm{~m})$。在直线轨迹与曲线轨迹的交界点，履带滑动参数均会突变。车辆按曲线路径行驶时，x_c增大，滑动参数R=(1.85 m, -1.63 m, 0.22 m), 这说明滑动参数与车辆状态和地形参数相关。

轨迹预测结果如图 7a所示，黑线代表车辆实际的行驶轨迹。图 7b与7c分别是图 7a中的局部放大图，轨迹预测误差如图 8所示。从图 7b与7c可以看出，EKF与UKF的轨迹预测误差相近，分别为0.18与0.16 m。3组不同预测时域参数的DUKF轨迹预测误差明显小于EKF与UKF。从图 7b与7c的局部放大图可以看出，预测时域参数为0.5 s时，DUKF的轨迹预测误差最小，相对应的预测轨迹更贴合参考轨迹，如图 8中的蓝线所示，其误差为0.03 m。预测时域参数为0.1与1 s的轨迹预测误差分别为0.10和0.06 m。对于小曲率工况，从轨迹预测误差与效果图中可以看出，过小的预测时域参数对上层UKF补偿较小，过大的预测时域参数会过多补偿上层UKF。选择合适的预测时域参数可以很大程度地提高参数估计的精度。

第2组参考轨迹设置为3段直线和2段曲率较大的圆弧，为大曲率工况。路径起点的坐标为(0, 0)，其曲率是连续的，如图 9所示。履带车辆两侧履带卷绕速度如图 10所示，转弯处的高速侧履带卷绕速度达到4.8 m/s，两侧履带卷绕速度差达0.7 m/s。

大曲率工况下，DUKF的滑动参数估计结果如图 11所示。在车辆刚起步时，y_l与y_r均处于两侧履带中心，在第1个弯道曲率最大处，滑动参数R=(2.72 m, -1.75 m, 0.25 m)。在第2个弯道曲率最大处，滑动参数R=(1.76 m, -2.72 m, 0.26 m)。

大曲率工况轨迹预测结果如图 12a所示，黑线代表车辆实际的行驶轨迹。轨迹预测误差如图 13所示。在2个弯道路径曲率最大处，UKF与EKF的轨迹预测误差分别大约为0.18、0.23 m(第1个弯道)与0.18、0.20 m(第2个弯道)。3组不同预测时域参数的DUKF轨迹预测误差明显小于EKF与UKF。图 12b为图 12a中曲率较大处的局部放大图，可以看出预测时域参数为0.5 s的预测轨迹更贴合参考轨迹，对应于图 13中52 s时的轨迹预测误差，为0.03 m左右。图 12c为图 12a中曲率较小处的局部放大图，预测时域参数为1 s的预测轨迹更贴合参考轨迹，对应于图 13中27 s时的轨迹预测误差，为0.04 m左右。在20与52 s处，下层UKF可以根据履带车辆的当前状态预测滑动参数的变化趋势，显著提高了算法在极限工况下的收敛速度与精度。

综上所述，小曲率工况的仿真结果与大曲率工况相似。从仿真结果的分析来看，对于下层UKF不同的预测时域参数，过小的预测时域参数对上层UKF补偿较小，过大的预测时域参数会过多补偿下层UKF。但本文所提出的DUKF在不同预测时域参数下仿真得到的轨迹预测误差均明显小于EKF和UKF，即DUKF滑动参数估计精度优于EKF和UKF。

4 结论

本文在常规UKF的基础上进一步提出了一种基于轨迹预测补偿的DUKF滑动参数实时估计方法，并通过MATLAB/Simulink与RecurDyn联合仿真验证了DUKF算法的有效性与实用性。该方法仅利用车辆传感器信息就能在路面先验条件未知的情况下估计履带滑动参数。在相同的工况下，与UKF和EKF的仿真结果进行对比，所提出的DUKF算法在大曲率与小曲率工况下均能大幅提升滑动参数估计精度。除此之外，在下层UKF中，对于大曲率路径与小曲率路径应该设置不同的预测时域参数来对上层UKF初步估计的滑动参数进行补偿，以进一步提升滑动参数估计的精度。