重载车辆轮胎模型参数辨识与灵敏度分析

黄通，刘志浩，高钦和，王冬，马栋

（火箭军工程大学导弹工程学院，陕西西安 710025）

重载子午轮胎作为重载车辆与地面直接接触的部件，其力学性能直接影响车辆动力学和行驶平顺特性.与常见轿车轮胎相比，重载子午轮胎具有气压高、阻尼低、花纹粗大、扁平率大的特点.研究分析重载子午轮胎的力学性能对设计和分析先进的底盘系统和优化先进车辆系统结构至关重要［1］.为了研究轮胎力学性能，国内外研究者提出了多种轮胎模型［2-3］.

魔术公式（Magic Formula，MF）轮胎模型以试验数据为基础，通过精确拟合试验数据，描述轮胎的力学性能.基于MF 模型，相关研究者已发展了许多适应于各种工况的不同改进形MF 模型［4-5］.PAC2002模型是Pacejka［6］不断更新和完善的一种广泛应用于车辆动力学仿真和分析的MF模型，能够研究不同载荷、倾角和胎压下的轮胎力学性能，具有较强的工程应用背景.但特征参数繁多，而且高度非线性，给特征参数的辨识带来很大困难.

目前，轮胎模型的特征参数辨识基本采用数值优化算法和智能搜索算法.遗传算法因为具有较强的鲁棒特性，最先开始被应用于轮胎模型的特征参数辨识［7-8］.遗传算法虽然能够在全局范围内逼近最优解，但存在局部搜索能力较差，收敛速度慢的问题.为此，相关研究者将全局范围逼近较强的智能搜索算法和局部搜索较强的数值优化算法结合起来，产生了多种轮胎参数辨识的混合优化算法，这些混合算法普遍提升了轮胎模型的辨识精度［9-11］.

基于此，本文选择适合在动态和多目标优化环境中寻优的粒子群算法和Levenberg-Marquardt（LM）算法进行纯纵滑和纯侧偏工况下的轮胎模型的特征参数辨识，并对辨识后的轮胎模型进行残差分析和结果验证，最后基于Sobol 灵敏度分析对该型轮胎进行参数分析和优化，研究结果可以为该型轮胎的应用提供理论支撑.

1 轮胎六分力测试试验

根据轮胎动力学原理，结合现有的试验设备，进行不同垂向载荷作用下纯纵滑和纯侧偏工况的轮胎六分力测试试验.试验轮胎型号为GL073A，试验平台为六自由度轮胎力学特性试验台，试验路面为4 m水泥路面滑台，滑台运行速度为200 mm∕s.纯纵滑工况试验中滑移率为-0.018～0；在纯侧偏工况试验中，侧偏角为-10°～13°.轮胎六分力测试试验如图1所示.

图1 轮胎六分力测试试验Fig.1 Six component force test for tire

在轮胎压力为810 kPa 下，当垂向载荷分别为24 800 N、34 700 N 和44 600 N 时，分别测试轮胎力学性能.测试试验结果分别如图2～图4所示.

图2 侧偏力与侧偏角关系曲线Fig.2 Curve of sideforce and sideslip angle

图3 回正力矩与侧偏角关系曲线Fig.3 Curve of aligning torque and sideslip angle

图4 纵向力与滑移率关系曲线Fig.4 Curve of longitudinal force and slip ratio

根据轮胎力学，侧偏刚度是侧偏角度为0 时的侧偏力曲线斜率；回正刚度是侧偏角度为0 时的回正力矩曲线斜率；纵滑刚度是滑移率为0 时的纵向力曲线斜率.由图2～图4可知：

1）当侧偏角不超过6°时，侧偏力与侧偏角呈线性关系，侧偏刚度从3 828 N∕（°）增加到5 320 N∕（°）；

2）侧偏角为6°时，回正力矩达到最大值，回正刚度从151 N·m∕（°）增加到287 N·m∕（°）；

3）纵滑刚度随垂向载荷增大，从301 799 增加到375 493.

2 魔术公式轮胎模型

魔术公式轮胎模型是基于试验数据通过一个公式即可表示轮胎的纵向力、侧偏力和回正力矩，公式中的特征参数物理意义明确.

基本的魔术公式为纯侧偏和纯纵滑工况下的侧偏力和纵向力：

式中：y为纯纵滑工况的纵向力或纯侧偏工况的侧偏力；x为对应的滑移率或侧偏角；B为刚度因子；C为形状因子；D为峰值因子；E为曲率因子.

对基本魔术公式进行修正，可以得到纯侧偏工况下的侧偏力计算公式：

式中：γ为外倾角；Fwz为轮胎垂向载荷；Fwz0为轮胎名义垂向载荷；dfwz为名义垂向载荷增量；Cy、Dy、Ey、Ky、By为侧偏工况侧偏力计算因子；SHy、SVy分别为侧偏工况横向和垂向漂移量；α为侧偏角；αy为修正侧偏角；Fy0为侧偏力；py为侧偏力特征参数.

在纯侧偏工况中，侧偏力计算公式中待辨识的特征参数共18个［12］，如表1所示.

表1 侧偏力的特征参数Tab.1 Characteristic parameters of sideforce

纯侧偏工况下的回正力矩可以用轮胎拖距和残余力矩来计算：

式中：t为轮胎拖距；Mzr为残余力矩.均可采用余弦函数表达.

式中：Ct、Dt、Et、Kt、Bt为侧偏工况回正力矩计算因子；SHt、SHr分别为侧偏工况拖距和残余力矩的横向漂移量；αt为拖距修正侧偏角；αr为残余力矩修正侧偏角；R0为轮胎半径；qz为特征参数.

式（6）中各因子均能由类似式（3）的计算公式表达出来，在纯侧偏工况中，回正力矩计算公式中待辨识的特征参数共25个［12］，如表2所示.

表2 回正力矩的特征参数Tab.2 Characteristic parameters of aligning torque

对基本魔术公式进行修正，可以得到纯纵滑工况下的纵向力计算公式：

式中：Cx、Dx、Ex、Kx、Bx为侧偏工况回正力矩计算因子；SHx、SVx分别为纵滑工况横向和垂向漂移量；λ为滑移率；λx为修正滑移率；Fx0为纵向力；px为特征参数.

在纯纵滑工况中，纵向力计算公式中待辨识的特征参数共15个［12］，如表3所示.

表3 纵向力的特征参数Tab.3 Characteristic parameters of longitudinal force

3 特征参数辨识

3.1 粒子群优化算法

粒子群优化算法的提出来源于鸟类群体行为的规律性，通过模拟鸟类的觅食行为，将寻求优解问题的搜索空间比作鸟类的飞行空间，每只鸟抽象成一个没有质量和体积的粒子，用该粒子表征问题的一个可能解，从搜索空间中的随机解出发，通过迭代运算寻找最优解.粒子群优化算法由于参数少且易于实现，在非线性、多峰问题中具有较强的全局搜索能力，在科学研究与工程实践中均得到广泛关注.

对于魔术公式轮胎模型的特征参数辨识问题，假设在一个D维特征参数向量解搜索空间中，有N个粒子组成一个群，其中第i个粒子在搜索空间中的位置可以表示为：

第i个粒子在搜索空间中的“飞行速度”可以表示为：

第i个粒子目前搜索到的最优解为：

整个粒子群目前搜索到的最优解为：

迭代过程中粒子更新后的粒子速度和位置分别为：

式中：v(t)、x(t)分别为目前粒子速度和位置；w为惯性权重；t为迭代次数；c1、c2为学习因子；r1、r2为［0，1］内的均匀随机数.

式（14）中惯性权重w表示粒子在多大程度保留原来的速度.相关试验研究发现，较大的惯性权重能够使算法全局收敛能力强，而较小的惯性权重能够使算法局部收敛能力强.因此，本文选择由Shi 等［13］提出的线性递减动态权值策略，其表达式如下：

式中：Tmax为最大进化代数；wmax、wmin分别为最大惯性权重和最小惯性权重.在大多数应用中，通常取wmax=0.9，wmin=0.4.

3.2 Levenberg-Marquardt优化算法

数值优化算法具有相对较强的局部搜索能力.本文将粒子群优化算法辨识的结果作为数值优化算法的初始值，选择收敛精度较高、非奇异性较好的Levenberg-Marquardt（LM）数值优化算法进行局部最优搜索.

LM 算法通常被视为Gauss-Newton 算法和最速下降法的结合，其基本原理公式为：

式中：x为特征参数向量；g(t)为残差；μ为参数因子；J为g(t)的Jacobian矩阵.

3.3 辨识结果

基于MATLAB 编程实现粒子群优化算法和LM优化算法.首先，利用粒子群算法的全局搜索能力辨识出魔术公式轮胎模型特征参数的近似解，并将其作为LM 算法的初始值进行局部精确搜索，得到特征参数的最终解.然后，采用式（18）对辨识结果进行残差定量分析与评价.

轮胎侧偏力、回正力矩和纵向力的辨识结果分别如图5～图7所示.

图5 侧偏力辨识结果Fig.5 Identification results of sideforce

图6 回正力矩辨识结果Fig.6 Identification results of aligning torque

图7 纵向力辨识结果Fig.7 Identification results of longitudinal force

由图5～图7 可以看出，通过粒子群优化算法的辨识结果与试验数据大致吻合，说明该算法能够有效获得魔术公式轮胎模型特征参数的近似解；通过混合优化算法的辨识结果与试验数据吻合度较高，辨识结果残差相对减小，说明混合优化算法能够获得更为精确的解.

受试验数据采集量多的影响，垂直载荷为44 100 N 时的辨识残差相对于垂直载荷分别为24 500 N 和34 300 N 时的辨识残差要低，辨识结果与试验数据的吻合度相对较高.由图7 可知，在纵向力辨识结果中，由于试验数据的类线性特征，混合优化算法的辨识结果对粒子群算法的辨识结果提升相对较弱.

4 多特征参数敏感性分析

轮胎模型特征参数灵敏度能够反映特征参数变化对轮胎模型响应的影响程度，灵敏度较大的特征参数对轮胎模型拥有更强的调整能力，通常可以当作主导特征参数以反映轮胎模型的变化情况，特别是魔术公式这类特征参数较多的轮胎模型，确定主导特征参数能够有效提高模型调整的便捷性.

Sobol 灵敏度分析方法是基于多重积分的分解方法，可以对非线性的复杂模型进行分析，研究各参数同时变化下各参数之间的耦合作用.本文选择Sobol 灵敏度分析方法对魔术公式轮胎模型特征参数的影响进行分析，并以各特征参数的一阶灵敏度和总阶灵敏度作为评价标准.

4.1 侧偏力特征参数

在采样数为5 000时，侧偏力特征参数的一阶和总阶灵敏度分析结果如图8所示.

图8 侧偏力特征参数灵敏度Fig.8 Sensitivity of characteristic parameters for sideforce

由图8 可知，垂直漂移因子特征参数pvy1的一阶灵敏度最大，约为0.35，总阶灵敏度稍大，约为0.18，说明pvy1不仅对模型响应影响较大，而且与其他参数存在耦合关系，这是因为在式（2）和式（3）中，pvy1作为一个叠加项独立存在，对模型影响比较直接.重点分析剩余一阶灵敏度较大的刚度因子特征参数pky1和pky2和水平漂移因子特征参数phy2.

在pky1和pky2以及phy2不同值下，各主要特征参数总阶灵敏度如图9所示.

图9 侧偏力主要特征参数的总阶灵敏度Fig.9 Distribution of total order sensitivity for sideforce main characteristic parameters

由图9可知：

1）峰值因子特征参数pdy1受pky1、pky2和phy2的影响较大，这与图8 所示的pdy1一阶灵敏度较小，总阶灵敏度较大的结果是一致的.其灵敏度随刚度因子特征参数pky1和pky2的增大而减小，随水平漂移因子特征参数phy2的增大而增大.

2）垂直漂移因子特征参数pvy1灵敏度随刚度因子特征参数pky1和pky2的增大而增大，随水平漂移因子特征参数phy2的增大而减小.

3）刚度因子特征参数pky1和pky2以及水平漂移因子特征参数phy2之间的耦合灵敏度较高，pky1和pky2相互耦合灵敏度呈正比关系.当phy2为14 时，pky1和pky2的灵敏度均达到最高.

在侧偏力18 个特征参数中，选择刚度因子特征参数pky1和pky2、水平漂移因子特征参数phy2和垂直漂移因子特征参数pvy1为主导特征参数.通过调整这4个参数来反映魔术公式轮胎模型侧偏力的变化规律，并选择峰值因子特征参数pdy1为主要特征参数，可以通过调整峰值因子特征参数，使主导特征参数更为准确.

4.2 回正力矩特征参数

在采样数为5 000时，回正力矩特征参数的一阶和总阶灵敏度分析结果如图10所示.

图10 回正力矩特征参数灵敏度Fig.10 Sensitivity of characteristic parameters for aligning torque

由图10 可知，一阶灵敏度大于0.10 的特征参数有qez1、qez4、qbz1和qhz1.其中，曲率因子特征参数qez1的一阶灵敏度最大，约为0.38，但总阶灵敏度相对较小，说明qez1对模型影响较大，但与其他参数耦合相对较小；曲率因子特征参数qez4和水平漂移因子特征参数qhz1也存在相同的规律.而峰值因子特征参数qdz1和形状因子特征参数qcz1虽然一阶灵敏度较小，但总阶灵敏度相对较大，qdz1的总阶灵敏度最大，约为3.75.因此，重点分析总阶灵敏度较大的qdz1和qcz1，如图11（a）和图11（b）所示.

图11 回正力矩主要特征参数的总阶灵敏度Fig.11 Distribution of total order sensitivity for aligning torque main characteristic parameters

由图11可知，峰值因子特征参数qdz1主要与形状因子特征参数qcz1耦合相关，二者相互耦合灵敏度较高，随着qdz1的增加，qcz1灵敏度增加幅度较小，基本控制在0.4～0.45，其他特征参数灵敏度也相对稳定.随着qcz1的增加，qdz1灵敏度增长幅度较大，当qcz1为-0.5 时，qdz1灵敏度达到最大，约为0.95；同时随着qcz1的增大，其他特征参数的灵敏度均逐渐减小.

在回正力矩的25 个特征参数中，选择曲率因子特征参数qez1和qez4、刚度因子特征参数qbz1、水平漂移因子特征参数qhz1以及峰值因子特征参数qdz1为主导特征参数；选择形状因子特征参数qcz1为主要特征参数.

4.3 纵向力特征参数

在采样数为5 000时，纵向力特征参数的一阶和总阶灵敏度分析结果如图12所示.

图12 纵向力特征参数灵敏度Fig.12 Sensitivity of characteristic parameters for longitudinal force

由图12 可知，刚度因子特征参数pkx1的一阶灵敏度和总阶灵敏度最大，分别约为0.37 和0.25，说明pkx1对模型响应和其他参数的耦合关系均有较大影响.剩余一阶灵敏度较大的有：水平漂移因子特征参数phx1、垂直漂移因子特征参数pvx1和峰值因子特征参数pdx1.其中pvx1为模型叠加项且总阶灵敏度较小，pdx1和phx1的总阶灵敏度都相对较大.因此，重点分析峰值因子特征参数pdx1和水平漂移因子特征参数phx1变化下的总阶灵敏度分布规律.

pdx1和phx1变化下的总阶灵敏度如图13 所示.由图13 可知，随着pdx1的增大，刚度因子特征参数pkx1和水平漂移因子特征参数phx1均表现出较高的灵敏度，并随之增加；其他特征参数灵敏度随着pdx1的增大而逐渐减小.随着水平漂移因子phx1的逐渐增大，pdx1和pkx1灵敏度虽然占据较大分量，但灵敏度总体呈减小趋势.

图13 纵向力主要特征参数的总阶灵敏度Fig.13 Distribution of total order sensitivity for longitudinal force main characteristic parameters

因此，在纵向力15 个特征参数中，选择刚度因子特征参数pkx1、水平漂移因子特征参数phx1、垂直漂移因子特征参数pvx1以及峰值因子特征参数pdx1为主导特征参数.

将通过Sobol 灵敏度分析选择出的主导特征参数代入轮胎模型进行计算，结果如表4所示.

表4 辨识结果对比Tab.4 Comparison of identification results

由表4 可知，采用Sobol 灵敏度分析所得的主导参数进行辨识的结果与直接采用混合优化进行辨识的结果相比，辨识残差均相对增大，但增大的幅值较小.其中垂向载荷为34 300 kN 的回正力矩残差增幅相对较大，残差增幅为0.138%，这是由于主导参数辨识省略了其他特征参数影响，使残差增大，但仍小于单独粒子群算法的辨识残差.而稳定迭代次数相对于直接采用混合优化进行辨识的结果减小明显，模型收敛速度增加，最大增幅为30.4%，对轮胎模型拥有较强的调整能力，采用Sobol 灵敏度分析所得出的主导参数能够有效提高模型调整的便捷性.

5 结论

本文以某型轮胎六分力测试试验为基础，对该型重载子午轮胎的魔术公式轮胎模型进行参数辨识和灵敏度分析.得出以下结论：

1）基于粒子群算法能够实现对魔术公式轮胎模型特征参数的辨识，且辨识结果残差分布在5%左右.混合优化算法能够提高特征参数辨识的精确性，将辨识结果残差控制在5%以内.

2）基于Sobol 灵敏度分析方法可以从魔术公式轮胎模型的58 个特征参数中选择出13 个特征参数作为主导特征参数，可以优先调整这13 个主导参数来控制魔术公式轮胎模型的变化规律.