重載車輛輪胎模型參數辨識與靈敏度分析

黃通，劉志浩，高欽和，王冬，馬棟

（火箭軍工程大學導彈工程學院，陜西西安 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 個主導參數來控制魔術公式輪胎模型的變化規律.