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

哈爾濱工業(yè)大學(xué) 機(jī)器人技術(shù)與系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150080

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

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

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

總體來(lái)說(shuō)，現(xiàn)有面向輪地力學(xué)模型的參數(shù)靈敏度分析以局部分析為主，且局限于定性層面，缺乏面向全局的定量分析研究。多數(shù)參數(shù)估計(jì)方法依賴于不同工況下的狀態(tài)數(shù)據(jù)，且以地面性質(zhì)均一穩(wěn)定為假設(shè)。不同工況通常需要人為變換車輪控制策略，難以滿足控制實(shí)時(shí)性需求，均質(zhì)地面假設(shè)在實(shí)際應(yīng)用中也較難滿足。因此，亟需面向輪地相互作用模型開展地面力學(xué)參數(shù)的全局靈敏度定量分析，篩選出主導(dǎo)參數(shù)并提出對(duì)應(yīng)的估計(jì)方法，為實(shí)時(shí)控制奠定基礎(chǔ)。

區(qū)別于定性的參數(shù)靈敏度分析和工況要求苛刻的全參數(shù)估計(jì)方法，本文的主要貢獻(xiàn)為：創(chuàng)新地引入了Sobol分析方法，定量地分析了輪地相互作用模型中地面承壓特性參數(shù)和剪切特性參數(shù)的全局靈敏度，為主導(dǎo)參數(shù)選取提供了更可靠的依據(jù)；推導(dǎo)了解析形式的主導(dǎo)參數(shù)估計(jì)模型，僅通過(guò)單一工況下的狀態(tài)數(shù)據(jù)估計(jì)主導(dǎo)參數(shù)以反映地面承壓和剪切特性變化。所提出的主參數(shù)估計(jì)方法適用于多種滑轉(zhuǎn)工況，且估計(jì)結(jié)果可用于牽引力預(yù)測(cè)，以提高星球車控制的安全性和穩(wěn)定性。

1 考慮輪刺效應(yīng)的輪地相互作用模型

為了增強(qiáng)星球車的移動(dòng)性能，通常在星球車車輪輪周安裝輪刺以提高其牽引能力。例如機(jī)遇號(hào)火星車，其車輪輪周均勻排布著截面較小的豎直型輪刺。星球車在平坦地面上平穩(wěn)運(yùn)動(dòng)時(shí)，輪地接觸區(qū)域的相互作用如圖1[12]所示。圖中，ω為車輪轉(zhuǎn)速;v為車輪前進(jìn)速度;z為車輪沉陷量;r為車輪半徑;h為車輪輪刺高度。

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

星球車的行駛速度較低，其穩(wěn)定行駛時(shí)各時(shí)刻車輪的運(yùn)動(dòng)狀態(tài)可視為準(zhǔn)靜態(tài)。不考慮車輪側(cè)偏，可通過(guò)應(yīng)力積分列寫平衡方程[12]：

(1)

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

rs=r+λsh

(2)

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

為了適應(yīng)車輪的形狀特征，以Bekker承壓模型[18]為基礎(chǔ)，正應(yīng)力可表示為輪地接觸角的函數(shù)：

(3)

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

(4)

θm=(c1+c2s)θ1

(5)

θ2=c3θ1

(6)

式中:s為滑轉(zhuǎn)率。

通常情況下θ2較小，可將其簡(jiǎn)化為0。

為了表現(xiàn)車輪滑轉(zhuǎn)沉陷的動(dòng)態(tài)效應(yīng)，沉陷指數(shù)N可表示為滑轉(zhuǎn)率s的線性函數(shù)[19]

N=n0+n1s

(7)

(8)

同時(shí)，考慮車輪輪刺對(duì)于剪切應(yīng)力的影響，以Janosi剪切模型[20]為基礎(chǔ)，土壤剪切應(yīng)力τ可表示為

(9)

式中：c、φ、K為土壤固有的剪切特性參數(shù)，分別表示內(nèi)聚力、內(nèi)摩擦角和剪切變形模量；j(θ)為土壤的剪切位移;σ(θ)為正應(yīng)力?？紤]輪刺效應(yīng)時(shí)的剪切位移[12]可表示為

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

(10)

(11)

(12)

式中：θ′1為考慮輪刺效應(yīng)的進(jìn)入角；Rj為考慮輪刺效應(yīng)的等效剪切半徑；sj1和sj2為過(guò)渡滑轉(zhuǎn)率。一般情況下，sj1=0.15，sj2=0.5。

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

2 地面力學(xué)特性參數(shù)靈敏度分析

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

2.1 Sobol靈敏度分析方法

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

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

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

(13)

式中：f0為常數(shù)項(xiàng)。其他各子項(xiàng)對(duì)其所包含每個(gè)元素變量的積分為0，即

(14)

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

(15)

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

(16)

(17)

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

(18)

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

(19)

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

(20)

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

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

(21)

(22)

(23)

(24)

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

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

(25)

(26)

(27)

本文分析各地面力學(xué)參數(shù)對(duì)輪地相互作用力和力矩的影響，使用各參數(shù)的一階靈敏度及其總階靈敏度作為評(píng)價(jià)指標(biāo)。

2.2 地面力學(xué)參數(shù)靈敏度分析

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

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

表1 地面力學(xué)參數(shù)取值區(qū)間Table 1 Value range of terrain mechanical parameters

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

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

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

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

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

3 主導(dǎo)參數(shù)估計(jì)

3.1 主導(dǎo)參數(shù)解析模型

在輪地相互作用模型中輪地作用力宏觀表現(xiàn)為土壤應(yīng)力的積分，在應(yīng)力公式的原始形式下，積分難以直接求解，因此需要對(duì)應(yīng)力的形式做適當(dāng)?shù)暮?jiǎn)化。參考常用的輪地相互作用模型的簡(jiǎn)化方法[25]，對(duì)沿接觸區(qū)域的輪地作用應(yīng)力分布進(jìn)行線性化處理

(28)

(29)

式中：σm為最大正應(yīng)力；τm為最大切應(yīng)力。

(30)

(31)

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

(32)

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

(33)

該假設(shè)在滑轉(zhuǎn)率較大的情況下合理成立。在大部分情況下，c1約為0.4，c2小于0.3[16]。又因?yàn)棣?約為0，θm接近θ1的一半。

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

(34)

在該主導(dǎo)參數(shù)解析模型中，沉陷指數(shù)N的解析表達(dá)式是關(guān)于{θ1,FN,MR,kc,kφ,r,b,rs}的函數(shù)，其中{θ1,FN,MR}是系統(tǒng)狀態(tài)參數(shù)，{kc,kφ}是非主導(dǎo)地面力學(xué)參數(shù)，{r,b,rs}是星球車的車輪尺寸參數(shù)；內(nèi)摩擦角φ的解析表達(dá)式則是關(guān)于{s,θ1,θ′1,FN,MR,c,K,r,b,rs}的函數(shù)，其中{s,θ1,θ′1,FN,MR}是系統(tǒng)狀態(tài)參數(shù)，{c,K}是非主導(dǎo)地面力學(xué)參數(shù)，{r,b,rs}是星球車的車輪尺寸參數(shù)。系統(tǒng)狀態(tài)參數(shù)可通過(guò)直接或間接測(cè)量得到，車輪尺寸參數(shù)與車輪設(shè)計(jì)相關(guān)可視為已知參數(shù)，而非主導(dǎo)參數(shù)則可根據(jù)地面力學(xué)特性賦予適當(dāng)?shù)牡湫椭担酝瓿蓪?duì)主導(dǎo)參數(shù)的估計(jì)。

3.2 基于解析模型的主導(dǎo)參數(shù)估計(jì)方法

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

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

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

算法1 主導(dǎo)參數(shù)估計(jì)輸入:沉陷量序列Sz={z1, z2, …, zn}, 滑轉(zhuǎn)率序列Ss={s1, s2, …, sn}, 法向力序列SF={FN1, FN2, …, FNn}, 驅(qū)動(dòng)阻力矩序列SM={MR1, MR2, …, MRn}, 元素?cái)?shù)量為n,非主導(dǎo)參數(shù){kc, kφ, c, K},車輪尺寸參數(shù){b, r, rs, h},濾波器寬度L,濾波器切換閾值{T1, T2, T3, T4},終止條件Ne輸出:沉陷指數(shù)序列SN, 內(nèi)摩擦角序列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 試驗(yàn)驗(yàn)證

4.1 基于仿真的地面力學(xué)特性主導(dǎo)參數(shù)估計(jì)試驗(yàn)

為了檢驗(yàn)?zāi)Ｐ偷臏?zhǔn)確性，通過(guò)仿真模擬輪地相互作用力和力矩進(jìn)行地面力學(xué)主導(dǎo)參數(shù)估計(jì)試驗(yàn)。仿真所用車輪相關(guān)尺寸參數(shù)參考星球車原理樣機(jī)的車輪[12]設(shè)置，各參數(shù)取值為：r=0.135 m，b=0.11 m，h=0.015 m。仿真模擬星球車車輪由地形I經(jīng)地形II再到地形I的直線行駛過(guò)程，地形I和地形II的地面力學(xué)參數(shù)[29]如表2所示。仿真所得系統(tǒng)狀態(tài)參數(shù)FN和MR通過(guò)添加高斯白噪聲以模擬實(shí)際測(cè)量值誤差，不同參數(shù)所添加的噪聲標(biāo)準(zhǔn)差列于表3中。各系統(tǒng)狀態(tài)參數(shù)的采樣頻率設(shè)為100 Hz。

表2 仿真地形的地面力學(xué)參數(shù)值Table 2 Mechanical parameters of terrain in simulation

表3 系統(tǒng)狀態(tài)參數(shù)噪聲標(biāo)準(zhǔn)差Table 3 Standard deviation of system state parameter noise

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

4.2 基于單輪試驗(yàn)的地面力學(xué)特性估計(jì)試驗(yàn)

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

表4 模擬月壤的地面力學(xué)參數(shù)值Table 4 Terrain mechanical parameters of lunar simulant

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

由于輪刺的作用，采集到的力和力矩?cái)?shù)據(jù)存在周期性的小幅波動(dòng)，根據(jù)對(duì)數(shù)據(jù)波形的分析，設(shè)置濾波器的寬度為3個(gè)波動(dòng)周期，濾波器方法切換的閾值為4倍波形幅值。利用采集的數(shù)據(jù)通過(guò)主導(dǎo)參數(shù)解析模型對(duì)沉陷指數(shù)N和內(nèi)摩擦角φ進(jìn)行估計(jì)。

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

4.3 滑轉(zhuǎn)率對(duì)估計(jì)結(jié)果的影響

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

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

4.4 非主導(dǎo)參數(shù)對(duì)估計(jì)結(jié)果的影響

在驗(yàn)證試驗(yàn)中，為了說(shuō)明主導(dǎo)參數(shù)解析模型的效果，采用了準(zhǔn)確的地面力學(xué)特性非主導(dǎo)參數(shù)值用于估計(jì)。而在實(shí)際應(yīng)用中，非主導(dǎo)參數(shù)通?；诮?jīng)驗(yàn)賦予標(biāo)稱值，該標(biāo)稱值可能與真實(shí)值之間存在偏差，因此需要分析非主導(dǎo)參數(shù)的變動(dòng)對(duì)于估計(jì)結(jié)果的影響。

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

由圖12可知，土壤摩擦變形模量kφ的減小，會(huì)使沉陷量指數(shù)N的估計(jì)值降低，并且所選用的kφ典型值與真實(shí)值偏差越大，對(duì)應(yīng)的沉陷指數(shù)估計(jì)值降低得越快，當(dāng)所選用的kφ典型值的與真實(shí)值的偏差在500 kPa/mN范圍內(nèi)時(shí)，N的估計(jì)誤差小于0.1。同時(shí)，隨著沉陷量的增加，kφ變動(dòng)對(duì)于N的影響稍有增加。對(duì)于內(nèi)摩擦角，隨著土壤剪切變形模量K的增大，φ的估計(jì)值也會(huì)相應(yīng)地增大，沉陷量為0.01 m時(shí)，K的值在0.003 m內(nèi)變化，φ的偏差小于4°；另外，隨著沉陷量的增加，內(nèi)摩擦角的估計(jì)結(jié)果受到土壤剪切變形模量的影響越小，K的影響就相應(yīng)減弱。當(dāng)沉陷量z=0.06 m 時(shí)，K的值增加一倍，其對(duì)于φ的影響小于6°。從總體來(lái)說(shuō)，土壤摩擦變形模量kφ和土壤剪切變形模量K對(duì)于主導(dǎo)參數(shù)估計(jì)結(jié)果的影響較小，并且當(dāng)取適當(dāng)?shù)牡湫椭禃r(shí)，估計(jì)結(jié)果可以保證有較高的精度。

5 結(jié) 論

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

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

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