輪地力學模型參數靈敏度分析與主參數估計

周如意，豐文浩，鄧宗全，高海波，丁亮，李楠

哈爾濱工業大學 機器人技術與系統國家重點實驗室，哈爾濱 150080

星球車作為典型的輪式移動機器人，不僅在深空探測領域發揮著關鍵作用，還是移動控制的重要研究對象。由于星球表面的地形松軟崎嶇，多輛星球車曾發生車輪打滑或沉陷事故。為了改善這種情況，提高星球車的移動性能，考慮輪地相互作用的控制算法[1-2]受到了廣泛的關注和研究。地面力學參數是輪地相互作用模型的重要組成部分，對星球車的控制效果具有直接影響。但由于載荷的限制，星球車幾乎無法搭載測量星壤力學特性的專用設備[3-4]，大部分地面力學參數無法通過直接測量獲得，因此地面力學參數的準確估計已成為限制控制算法應用的關鍵。主流的地面力學參數估計方法以輪地相互作用模型為基礎，但該模型形式復雜，參數眾多，通常無法同時對所有參數進行準確辨識。因此，對輪地相互作用影響較大的主導參數進行實時估計，以反映地形變化情況，對非主導參數進行在線調整，以提高參數精度，成為支持實時控制的一種折衷辦法。

為了篩選出模型的主導參數，通常需要分析參數的靈敏度。在地面力學領域，一般采用局部分析法求解地面力學參數的靈敏度。Li等[5]采用直接求導法分析了地面承壓特性模型和剪切特性模型中參數的靈敏度，結合參數取值獲得模型的主導參數。但局部分析法只能求解參數在值域內某點附近的局部靈敏度，當輸入變量的變化處于不同量級時或模型具有非線性特性時，對應的靈敏度分析能力有限。而全局靈敏度分析方法對分析較復雜的模型具有較大優勢，已有多種全局靈敏度分析方法被提出，如Morris方法[6]、Sobol方法[7]、傅里葉幅值敏感性測試[8]等，并廣泛應用于水文模型分析[9]和工程結構分析[10]等領域。鄧宗全等[11]采用數值法分析了車輪沉陷量，驅動阻力矩和掛鉤牽引力對土壤參數變化的敏感程度，得出土壤承壓特性參數主要對車輪沉陷量產生影響，土壤剪切特性參數對驅動阻力矩和掛鉤牽引力的影響較大，輪地作用對于接觸角系數的變化不敏感等定性的結論。

由于地面力學參數較多，利用輪地作用力學平衡方程一般無法直接求解所有參數值，通常利用不同工況(例如不同滑轉率)下的車輪狀態數據迭代優化求解。Ding等[12]推導了封閉形式的地面力學參數解析解耦方程，并提出將地面力學的8個參數分為3個階段逐步進行求解的在線辨識方法。Hutangkabodee等[13]使用復合辛普森法則代替輪地相互作用模型中應力的直接積分，簡化了輪地相互作用模型，并采用廣義牛頓-拉夫森方法實現了地形參數的實時辨識，所提出的方法在車輪牽引力預測中得到運用。Song等[14]根據車輪受到的垂直載荷和電機驅動力矩，采用遞歸神經網絡分別對地面承壓特性和剪切特性參數進行了自適應辨識。Xue等[15]使用了多輸出最小二乘支持向量機，在不需要測量車輪沉陷量的情況下，利用車輪滑轉率、載荷和轉矩等信息完成了地形剪切特性參數的在線預測。為了進一步提高參數估計的實時性，適應實時控制的需求，部分工作通過推導解析模型進行實時參數估計。Iagnemma等[16]通過線性化應力公式簡化了輪地作用模型，并推導出松軟地面下土壤內聚力和內摩擦角的解析表達式，用于在線計算土壤參數。Li等[5]利用三角函數擬合車輪的應力分布以簡化傳統的輪地相互作用模型，并提出了一種具有容錯開關的魯棒卡爾曼濾波器用于模型中的主要參數實時估計，所提出的估計方法能夠快速適應土壤類型的變化。其后在此基礎上又提出了地形參數估計的雙層結構[17]，內層實現主導參數的實時更新，外層采用遞歸高斯-牛頓算法對所有參數進行調整，提高了估計結果的精度。

總體來說，現有面向輪地力學模型的參數靈敏度分析以局部分析為主，且局限于定性層面，缺乏面向全局的定量分析研究。多數參數估計方法依賴于不同工況下的狀態數據，且以地面性質均一穩定為假設。不同工況通常需要人為變換車輪控制策略，難以滿足控制實時性需求，均質地面假設在實際應用中也較難滿足。因此，亟需面向輪地相互作用模型開展地面力學參數的全局靈敏度定量分析，篩選出主導參數并提出對應的估計方法，為實時控制奠定基礎。

區別于定性的參數靈敏度分析和工況要求苛刻的全參數估計方法，本文的主要貢獻為：創新地引入了Sobol分析方法，定量地分析了輪地相互作用模型中地面承壓特性參數和剪切特性參數的全局靈敏度，為主導參數選取提供了更可靠的依據；推導了解析形式的主導參數估計模型，僅通過單一工況下的狀態數據估計主導參數以反映地面承壓和剪切特性變化。所提出的主參數估計方法適用于多種滑轉工況，且估計結果可用于牽引力預測，以提高星球車控制的安全性和穩定性。

1 考慮輪刺效應的輪地相互作用模型

為了增強星球車的移動性能，通常在星球車車輪輪周安裝輪刺以提高其牽引能力。例如機遇號火星車，其車輪輪周均勻排布著截面較小的豎直型輪刺。星球車在平坦地面上平穩運動時，輪地接觸區域的相互作用如圖1[12]所示。圖中，ω為車輪轉速;v為車輪前進速度;z為車輪沉陷量;r為車輪半徑;h為車輪輪刺高度。

車輪所受到的外力包括豎直方向的垂直載荷W，水平方向的前進阻力fDP，電機的驅動力矩T和土壤對車輪的作用力。松軟土壤對車輪的作用力表現為連續的正應力σ和剪應力τ，分布情況受到車輪進入角θ1和車輪離去角θ2的影響，且在θm處產生最大應力。由于輪刺的存在，當車輪轉動時，車輪輪周表面與土壤的摩擦轉變為土壤與土壤之間的摩擦，并可等效為以等效剪切半徑rs(圖1中虛線所示)作用于車輪，對應的剪切應力相較于光輪情況顯著增加。應力的微觀作用通過積分表現為宏觀上的力和力矩，即法向力FN、掛鉤牽引力FDP和驅動阻力矩MR。

星球車的行駛速度較低，其穩定行駛時各時刻車輪的運動狀態可視為準靜態。不考慮車輪側偏，可通過應力積分列寫平衡方程[12]：

(1)

式中：b為車輪寬度；θ為輪地作用角；σ1和τ1分別為車輪前部分的正應力和切應力；σ2和τ2分別為車輪后部分的正應力和切應力。當輪刺高度為h時，等效的剪切半徑rs可表示為

rs=r+λsh

(2)

式中：λs為輪刺系數，通常取0.5[12]。

為了適應車輪的形狀特征，以Bekker承壓模型[18]為基礎，正應力可表示為輪地接觸角的函數：

(3)

式中：kc、kφ、N是描述土壤承壓特性的固有參數，分別為土壤內聚變形模量、摩擦變形模量和沉陷指數。輪地接觸角θ1、θ2和θm是沉陷量與應力分布系數系數c1、c2和c3的函數。

(4)

θm=(c1+c2s)θ1

(5)

θ2=c3θ1

(6)

式中:s為滑轉率。

通常情況下θ2較小，可將其簡化為0。

為了表現車輪滑轉沉陷的動態效應，沉陷指數N可表示為滑轉率s的線性函數[19]

N=n0+n1s

(7)

(8)

同時，考慮車輪輪刺對于剪切應力的影響，以Janosi剪切模型[20]為基礎，土壤剪切應力τ可表示為

(9)

式中：c、φ、K為土壤固有的剪切特性參數，分別表示內聚力、內摩擦角和剪切變形模量；j(θ)為土壤的剪切位移;σ(θ)為正應力。考慮輪刺效應時的剪切位移[12]可表示為

j(θ)=rs[(θ′1-θ)-(1-s)(sinθ′1-sinθ)]

(10)

(11)

(12)

式中：θ′1為考慮輪刺效應的進入角；Rj為考慮輪刺效應的等效剪切半徑；sj1和sj2為過渡滑轉率。一般情況下，sj1=0.15，sj2=0.5。

根據以上分析，輪地相互作用受到地面力學特性的影響，分別由法向承壓特性參數組{kc,kφ,N}和切向剪切特性參數組{c,φ,K}具體表征。通過地面力學參數利用輪地相互作用模型可計算輪地相互作用力和力矩；反之亦可利用輪地相互作用力和力矩進行地面力學特性參數的估計。

2 地面力學特性參數靈敏度分析

參數的靈敏度反映了參數變化對數學模型響應的影響程度，靈敏度較大的參數對模型擁有更強的調整能力，可視為主導參數用于直接反映模型的變化情況。對輪地相互作用模型中的參數進行靈敏度分析能夠確定模型的主導參數，在應用中可僅用主導參數估計結果反映相互作用力變化，以提高控制的實時性。輪地相互作用模型中含有積分公式，為非線性模型，參數間還存在較強的耦合性，用局部法分析各參數的靈敏度過程復雜且精度不高。因此，根據輪地相互作用模型的特點，采用Sobol方法對模型中的地面力學特性參數進行全局靈敏度分析。

2.1 Sobol靈敏度分析方法

Sobol靈敏度分析方法[21]是基于多重積分的分解方法，可以對非單調、非線性、非疊加等復雜模型進行分析，且允許分析各參數同時變化的情況，用于獲得參數之間的耦合作用，參數的變化范圍可拓展到整個參數定義域。

Sobol靈敏度分析方法假設模型中需要分析的參數個數為n，根據各參數的取值范圍，定義一個n維單元體In={x|0≤xi≤1;i=1,2,…,n}作為輸入參數的空間域。Sobol方法的核心思想是把平方可積函數f(x)分解為如式(13)所示的常數項、單個參數以及各參數間相互結合的函數項[22]

+f12…n(x1,x2,…,xn)

(13)

式中：f0為常數項。其他各子項對其所包含每個元素變量的積分為0，即

(14)

將式(13)的左右兩邊分別在整個參數空間域內平方并積分，結合式(14)可以得到

(15)

則函數f(x)的總方差可以表示為

(16)

(17)

根據式(17)函數總方差等于各階方差之和，有

(18)

f(x)的各階偏方差與其總方差的比值就是參數的各階靈敏度，則全局靈敏度指數Si1i2…is為

(19)

同時，根據式(19)可得所有參數的靈敏度之和為1，即

(20)

在Sobol方法中，低階靈敏度反映了單個參數對模型輸出的影響，高階靈敏度反映了參數間的耦合關系。由式(19)推導出的Si是參數xi的一階靈敏度，表示參數xi發生變化時對函數f(x)的影響。Sij(i≠j)為二階靈敏度，反映參數xi與xj同時變化時對輸出的影響，同時考慮2個參數之間的耦合關系。

根據以上的理論分析，在實踐中利用蒙特卡洛法進行數值積分代替解析積分估計輸出結果的方差，當模型參數樣本數足夠大時，計算結果逼近解析解。因此可由以下公式進行靈敏度分析

(21)

(22)

(23)

(24)

式中：xim為求解域In空間中的采樣點；l為蒙特卡洛采樣數;上標(1)和(2)表示參數xim的2個l×n維求解域內的采樣數組;x(～i)m表示除參數xim外的參數，即

x(～i)m=(x1m,…,x(i-1)m,x(i+1)m,…,xnm)

(25)

(26)

(27)

本文分析各地面力學參數對輪地相互作用力和力矩的影響，使用各參數的一階靈敏度及其總階靈敏度作為評價指標。

2.2 地面力學參數靈敏度分析

采用Sobol全局靈敏度分析方法，分別以輪地相互作用力和力矩為研究對象，對地面力學參數進行靈敏度分析。輪地作用模型中車輪所受土壤的宏觀作用力在微觀表現為正應力和切應力的組合積分，因此可分別從表征地面正應力和切應力的力學參數中篩選出主導參數，即分別對承壓特性參數kc，kφ和N，剪切特性參數c，φ和K進行靈敏度比較。通常情況下，應力最大值位于進入角和離去角中間，應力分布系數c1和c2變動范圍不大且受到車輪形狀的影響，因此不對其做靈敏度分析，一般在研究中給定典型值。本文在試驗過程中，應力分布系數的取值分別為c1=0.5，c2=0。

對地面力學參數利用Sobol方法進行靈敏度分析需要確定各參數的定義域。根據多種土壤各地面力學參數的測量值[13,16,24-25]確定取值范圍，如表1所示。

表1 地面力學參數取值區間Table 1 Value range of terrain mechanical parameters

輪地相互作用模型中的其他參數設置如下:車輪寬度b=0.15 m，車輪半徑r=0.14 m，等效剪切半徑rs=0.145 m，滑轉率s=0.3。

采用Sobol方法對參數進行靈敏度分析需要選擇合適的采樣數，采樣數過少時靈敏度分析結果無法收斂，采樣數過大則浪費計算資源。這里通過分析沉陷量z=0.03 m時各參數對掛鉤牽引力的靈敏度隨樣本數的變化情況用于確定合適的采樣數。對承壓特性參數進行分析時，控制剪切特性參數不變，參數取值分別為c=2 kPa，φ=30°，K=0.015 m。同理，在進行剪切特性參數分析時，保持承壓特性參數不變，參數取值分別為kc=5 kPa/mN-1，kφ=1 500 kPa/mN，N=1。

不同采樣大小下的各地面力學參數的一階及總階靈敏度分析結果如圖2中所示。由曲線可以看出，在樣本數小于500時各參數的靈敏度存在較大的波動，樣本數大于1 000后靈敏度分析結果趨于穩定。因此，在后續試驗中設置樣本數為1 500，可使得靈敏度分析結果收斂并得出具有代表性的結論。整體變化趨勢表明樣本數對參數的靈敏度影響不大。其中，沉陷指數N始終是承壓特性中靈敏度最大的參數，其一階靈敏度約為0.7，總階靈敏度稍大，約為0.9，說明沉陷指數與其他參數存在耦合關系；在剪切特性參數中，內摩擦角靈敏度最大，一階靈敏度和總階靈敏度約為0.8。

地面承壓模型中正應力為沉陷量的函數，剪切模型中切應力又與正應力有關，因此參數靈敏度的分析需要考慮沉陷量的影響。使車輪的沉陷量在0.01～0.06 m之間變化，利用輪地相互作用模型計算各參數的靈敏度。根據前面的分析，設置樣本數量為1 500。不同沉陷量下，各地面力學參數對于車輪受到的法向力、掛鉤牽引力和驅動阻力矩的總階靈敏度如直方圖3(a)～圖3(c)所示。承壓特性參數方面，總體上各地面承壓特性參數的靈敏度受沉陷量影響較小，沉陷指數N相對于法向力FN、掛鉤牽引力FDP和驅動阻力矩MR的靈敏度均大于0.8。剪切特性參數方面，在沉陷量為0.01 m時，土壤剪切變形模量K和內摩擦角φ都表現出較高的靈敏度，內摩擦角的靈敏度略大。隨著沉陷量的增加，內摩擦角的靈敏度逐漸增大，土壤剪切變形模量的靈敏度相對減小，在沉陷量為0.02 m時，內摩擦角相對于各力和力矩的靈敏度約增大至剪切變形模量的2倍。

綜合以上靈敏度分析結果，可以得出：① 在表征地面承壓特性的參數中，沉陷指數的靈敏度最高，總階靈敏度高于0.8，其對正應力的影響最大；② 在表征地面剪切特性的參數中，內摩擦角對于切應力的影響在總體上來說最大，同時土壤剪切變形模量在沉陷量較小時對于切應力的影響也不可忽略。因此，選擇沉陷指數N和內摩擦角φ作為輪地相互作用模型的主導參數，通過估計這2個參數反映地面在承壓和剪切特性方面的變化情況；其他參數(土壤內聚變形模量、土壤摩擦變形模量、土壤內聚力和土壤剪切變形模量)作為非主導參數，可根據經驗賦予典型值，以提供更準確的主導參數估計結果。

3 主導參數估計

3.1 主導參數解析模型

在輪地相互作用模型中輪地作用力宏觀表現為土壤應力的積分，在應力公式的原始形式下，積分難以直接求解，因此需要對應力的形式做適當的簡化。參考常用的輪地相互作用模型的簡化方法[25]，對沿接觸區域的輪地作用應力分布進行線性化處理

(28)

(29)

式中：σm為最大正應力；τm為最大切應力。

(30)

(31)

jm=rs[θ′1-θm-(1-s)(sinθ′1-sinθm)]

(32)

式中：jm為最大剪切位移。采用該線性化方法計算應力的簡化誤差小于15%[26]。其中最大應力角θm可表示為

(33)

該假設在滑轉率較大的情況下合理成立。在大部分情況下，c1約為0.4，c2小于0.3[16]。又因為θ2約為0，θm接近θ1的一半。

利用輪地相互作用平衡方程估計主導參數，由于待辨識的主導參數僅有2個，而輪地相互作用平衡方程有3個，存在冗余，可選擇其中2個方程用于主導參數求解。在實際應用中，為了減輕載荷，星球車車輪上一般不配備多維力傳感器，無法直接獲得車輪所受力與力矩數據。星球車的正常行駛速度較低，車輪的垂直載荷W可通過準靜態分析求解得到，車輪的電機驅動力矩T可由車輪電機的輸入電流估計得到，而前進阻力fDP不易直接獲取或求解。因此，選擇車輪垂直載荷W和電機驅動力矩T的平衡方程(式(1))用于參數求解，將簡化后的正壓力和切應力代入其中，經過積分與變換合并，可推導出主導參數解析表達式為

(34)

在該主導參數解析模型中，沉陷指數N的解析表達式是關于{θ1,FN,MR,kc,kφ,r,b,rs}的函數，其中{θ1,FN,MR}是系統狀態參數，{kc,kφ}是非主導地面力學參數，{r,b,rs}是星球車的車輪尺寸參數；內摩擦角φ的解析表達式則是關于{s,θ1,θ′1,FN,MR,c,K,r,b,rs}的函數，其中{s,θ1,θ′1,FN,MR}是系統狀態參數，{c,K}是非主導地面力學參數，{r,b,rs}是星球車的車輪尺寸參數。系統狀態參數可通過直接或間接測量得到，車輪尺寸參數與車輪設計相關可視為已知參數，而非主導參數則可根據地面力學特性賦予適當的典型值，以完成對主導參數的估計。

3.2 基于解析模型的主導參數估計方法

基于提出的主導參數解析模型，固定非主導參數，可對主導參數進行估計。該方法本質上是將主導參數視為模型整體性質的外在表現，主導參數的變動情況可以直接反映土壤性質的變動，從而使星球車可以根據土壤性質的改變調整其控制和規劃策略以提高其移動性能。具體估計方法如圖4所示。

車輪所受的法向力FN和驅動阻力矩MR通過車輪垂直載荷W和電機驅動力矩T反映，且分別通過車體準靜態分析和電機電流估計得到。車輪沉陷量z可通過外部視覺系統[27]或車體的位姿進行估計，滑轉率s通過測量車輪轉動角速度和行駛的線速度求解得到。車輪的角速度則通過車輪中安裝的編碼器測量，線速度可通過慣性導航原件測量[26-27]或通過視覺里程計[28]和車體運動學求解。星球車車輪尺寸參數已知，中間參數θ1，θ′1和rs可由r，z，h間接計算得到。另外，非主導地面力學參數可根據地形賦予典型估計值。

由于系統狀態參數的測量值存在噪聲，因此對數據進行濾波處理以減弱噪聲對估計結果準確性的影響。考慮星球車不同的運行狀態，當其穩定低速行駛時采用遞推平均濾波方法；當車輪經過不同地形交界處時，由于地面的力學特性發生變化，對應的系統狀態參數將隨之劇烈變化，若變化幅度超過設定閾值，則采用中值濾波對原始數據進行處理。采用均值濾波可以抑制車體平穩運行時的數據噪聲，同時起到平滑的作用；在地形變化處利用中值濾波的方法能夠保護數據尖銳的邊緣信號，實現數據的快速改變；均值濾波與中值濾波相比計算復雜度低，可在車體平穩運行時加快計算速度。對應偽代碼如算法1所示。

算法1 主導參數估計輸入:沉陷量序列Sz={z1, z2, …, zn}, 滑轉率序列Ss={s1, s2, …, sn}, 法向力序列SF={FN1, FN2, …, FNn}, 驅動阻力矩序列SM={MR1, MR2, …, MRn}, 元素數量為n,非主導參數{kc, kφ, c, K},車輪尺寸參數{b, r, rs, h},濾波器寬度L,濾波器切換閾值{T1, T2, T3, T4},終止條件Ne輸出:沉陷指數序列SN, 內摩擦角序列Sφ1: fori←1 tondo2: ^zi←DATAFILTER({z1, z2, …, zi}, L, T1)3: ^si←DATAFILTER({s1, s2, …, si}, L, T2)4: ^FNi←DATAFILTER({FN1, FN2, …, FNi}, L, T3)5: ^MRi←DATAFILTER({MR1, MR2, …, MRi}, L, T4)6: θ1←θ1(^zi, r)7: θ'1←θ'1(^zi, ^si, r, h)8: ^Ni←N(θ1, ^FNi, ^MRi, kc, kφ, r, b, rs)9: ^φi←φ(^si,θ1,θ'1,^FNi, ^MRi, c, K, r, b, rs)10: i←i + 111: add ^Ni to SN

12: add ^φi to Sφ13:end for14:15:Function DATAFILTER(S, L, T)16: i←number of elements in sequence S17: ai←i-th element in sequence S18: ai-1←(i-1)-th element in sequence S19: flag L20: ifiTthen24: flag←025: end if26: if flag < L then27: afiltered←the median of sequence {ai-L+1, ai-L+2,…,ai}28: flag←flag+130: else31: afiltered←the mean of sequence {ai-L+1, ai-L+2,…,ai}33: end if34: returnafiltered34: end if35:end function

4 試驗驗證

4.1 基于仿真的地面力學特性主導參數估計試驗

為了檢驗模型的準確性，通過仿真模擬輪地相互作用力和力矩進行地面力學主導參數估計試驗。仿真所用車輪相關尺寸參數參考星球車原理樣機的車輪[12]設置，各參數取值為：r=0.135 m，b=0.11 m，h=0.015 m。仿真模擬星球車車輪由地形I經地形II再到地形I的直線行駛過程，地形I和地形II的地面力學參數[29]如表2所示。仿真所得系統狀態參數FN和MR通過添加高斯白噪聲以模擬實際測量值誤差，不同參數所添加的噪聲標準差列于表3中。各系統狀態參數的采樣頻率設為100 Hz。

表2 仿真地形的地面力學參數值Table 2 Mechanical parameters of terrain in simulation

表3 系統狀態參數噪聲標準差Table 3 Standard deviation of system state parameter noise

對仿真數據進行主導參數估計，所用濾波器的寬度L設置為11，濾波方法切換的閾值Thx為5倍白噪聲的標準差。經過濾波處理和未經濾波處理的主導參數估計結果對比如圖5所示。從圖中可以看出，對系統狀態參數測量值進行濾波可明顯改善噪聲對主導參數估計結果的影響，在第一種地形中沉陷指數和內摩擦角估計結果的標準差分別由1.97×10-3和0.117°減小為6.20×10-4和0.036 3°。同時，所設計的濾波方法能夠快速跟隨地形的變化，對地形交接處突變數據的估計結果也能保持較好的波形。對地形I，沉陷指數的平均估計誤差為3.11%，內摩擦角的平均估計誤差為0.03%。對于地形II，相應的誤差分別是3.12%和0.11%。在同一地形中采用遞推平均濾波代替中值濾波有效地降低了計算量。因此，主導參數估計方法能夠快速跟隨地面特性的變化，同時較大程度上抑制測量噪聲，獲得較為準確的估計結果。

4.2 基于單輪試驗的地面力學特性估計試驗

為了檢驗主導參數解析模型及其估計方法在真實輪地相互作用中的效果，使用哈爾濱工業大學研制的星球車單輪測試平臺(如圖6所示)進行不同工況下的輪地相互作用試驗進行驗證。使用單個星球車車輪在平坦的模擬月壤[29]中進行輪地相互作用試驗，所用模擬月壤經過通風干燥、過篩和烘烤等處理，其地面力學參數已在表4中列出。所用車輪半徑r為0.157 3 m，車輪寬度b為0.165 m，車輪輪周均布24個高度h為0.015 m的直型輪刺。

表4 模擬月壤的地面力學參數值Table 4 Terrain mechanical parameters of lunar simulant

單輪測試臺驅動車輪運動，車輪的實際速度由拖拽電機控制，角速度由驅動電機控制，通過改變2個電機的轉速實現不同滑轉率的工況。車輪上方安裝有直線位移傳感器用于測量車輪沉陷量。車輪與支架連接處安裝有六維力傳感器和扭矩傳感器，能夠測量車輪垂直載荷、前進阻力及電機驅動力矩。實驗中，車輪輪刺的存在使得載荷波動范圍較大，為了減小載荷波動帶來的相對誤差，將垂直載荷設定為比較有代表性的80 N。參考實際星球車的移動速度，設定車輪的線速度為0.01 m/s，為避免車輪發生過大的沉陷，設定車輪滑轉率在0～0.6之間變化。驅動車輪在模擬月壤I中進行直線運動實驗，在車輪達到穩定狀態后采集力、力矩和沉陷量等車輪狀態數據，系統的采樣頻率為6.67 Hz。具體試驗過程如圖7所示。

由于輪刺的作用，采集到的力和力矩數據存在周期性的小幅波動，根據對數據波形的分析，設置濾波器的寬度為3個波動周期，濾波器方法切換的閾值為4倍波形幅值。利用采集的數據通過主導參數解析模型對沉陷指數N和內摩擦角φ進行估計。

滑轉率為0.4時的地面力學主導參數估計結果如圖8所示。在輪刺造成波動的前3個周期內，測量得到的參數不滿足濾波條件沒有進行濾波，故估計結果存在較大的波動；第3個周期之后，沉陷指數和內摩擦角的估計結果分別約在1.17和32°附近小幅波動。估計得到的地面力學主導參數以及系統狀態參數可代入式(1)中預測車輪受到的掛鉤牽引力，預測值與實際測量值如圖9所示。測量值在17～25 N間波動，預測值約在20 N附近，表明利用求解得到的地面力學主導參數能夠較準確地預測地面給車輪的牽引力。

4.3 滑轉率對估計結果的影響

此外，本文還利用不同滑轉率下車輪在模擬月壤II中的輪地相互作用實驗檢驗了該估計方法的滑轉率適用范圍。車輪滑轉率在0～0.6之間的地面力學特性主導參數估計結果如圖10所示。從圖中可以得出，沉陷指數與滑轉率呈線性關系，這也說明了將沉陷指數表示為N=n0+n1s的形式的合理性。沉陷指數預測值與真實值的誤差非常小，不同滑轉率下的平均相對誤差為2.8%。內摩擦角在滑轉率0～0.6之間的估計值與真實值偏差較小，平均相對誤差為2.8%。

利用主導參數的估計值預測車輪的掛鉤牽引力(圖11)，在滑轉率小于0.3時，預測值與測量值相差較小，隨著滑轉率增大，兩者之間的偏差逐漸增大，平均滿量程誤差為7.4%。在主導參數估計和掛鉤牽引力預測過程中，測量值與預測值之間較大的偏差主要由固化輪地接觸角系數c1，c2和c3的參數值引起，在僅以有限主導參數擬合模型的情況下，對掛鉤牽引力7.4%的平均滿量程誤差結果在可接受范圍內。

4.4 非主導參數對估計結果的影響

在驗證試驗中，為了說明主導參數解析模型的效果，采用了準確的地面力學特性非主導參數值用于估計。而在實際應用中，非主導參數通常基于經驗賦予標稱值，該標稱值可能與真實值之間存在偏差，因此需要分析非主導參數的變動對于估計結果的影響。

根據前文地面承壓特性參數和剪切特性參數的靈敏度分析結果，得出靈敏度僅次于沉陷指數N和內摩擦角φ的參數分別為土壤摩擦變形模量kφ和土壤剪切變形模量K，而土壤內聚變形模量kc和土壤內聚力c的靈敏度都接近于0，其對于主導參數求解結果的影響可以忽略。為了驗證估計方法對于非主導參數變化的魯棒性，變換不同的kφ和K進行估計試驗。同時為了說明車輪沉陷量的影響，控制沉陷量在0.01～0.06 m內改變，估計結果如圖12所示。

由圖12可知，土壤摩擦變形模量kφ的減小，會使沉陷量指數N的估計值降低，并且所選用的kφ典型值與真實值偏差越大，對應的沉陷指數估計值降低得越快，當所選用的kφ典型值的與真實值的偏差在500 kPa/mN范圍內時，N的估計誤差小于0.1。同時，隨著沉陷量的增加，kφ變動對于N的影響稍有增加。對于內摩擦角，隨著土壤剪切變形模量K的增大，φ的估計值也會相應地增大，沉陷量為0.01 m時，K的值在0.003 m內變化，φ的偏差小于4°；另外，隨著沉陷量的增加，內摩擦角的估計結果受到土壤剪切變形模量的影響越小，K的影響就相應減弱。當沉陷量z=0.06 m 時，K的值增加一倍，其對于φ的影響小于6°。從總體來說，土壤摩擦變形模量kφ和土壤剪切變形模量K對于主導參數估計結果的影響較小，并且當取適當的典型值時，估計結果可以保證有較高的精度。

5 結 論

1) 在表征地面承壓特性的參數中，沉陷指數N的總階靈敏度高于0.8，在地面承壓特性參數中靈敏度最高，對輪地作用力和力矩的影響最大；內摩擦角對于輪地作用力和力矩的影響在地面剪切特性參數中最大。因此，沉陷指數N和內摩擦角φ可作為主導參數用于反映地面承壓和剪切特性的變化。

2) 所推導的地面力學主導參數的解析模型以輪地相互作用力和力矩為輸入，以地面力學特性主導參數為輸出，在已知車輪尺寸參數的情況下，以標稱值固定非主導參數的方式，可用于地面力學特性主導參數估計。

3) 所提出的主導參數估計方法通過對系統狀態參數進行濾波以抑制測量噪聲能夠快速跟隨地形的變化。在不同滑轉率下，主導參數估計結果總體上表現出較小的誤差，沉陷指數估計值較真實值的平均相對誤差為2.8%，內摩擦角的平均相對誤差小于3%，且估計結果對非主導參數的變化具有魯棒效果。本文的估計方法對精確的全參數估計和基于輪地相互作用的控制均具有參考價值。