王志福， 王學晨， 王 陽， 梁常春， 王 瑞

（1.北京理工大學電動車輛國家工程研究中心，北京 100081；2.空間智能機器人系統(tǒng)技術與應用北京市重點實驗室，北京空間飛行器總體設計部，北京 100094）

0 引言

載人月球探測技術既是無人探月技術發(fā)展的延續(xù)，又是建造月球基地目標實現(xiàn)的基礎，標志著月球探測開始進入長期發(fā)展階段［1-2］。作為載人月球探測中不可或缺的角色，載人月球車是宇航員能在月面順利完成各項探測任務的保證。月球表面地形信息對載人月球車而言是十分重要的，尤其是獲取準確的月壤參數(shù)信息在很大程度上影響著月球車的行駛性能。對月壤參數(shù)的準確估計，有利于后續(xù)對月球車行駛過程中運動狀態(tài)量進行控制，進一步優(yōu)化月球車行駛能力，提高月球車操縱穩(wěn)定性，也能為未來探月技術研究發(fā)展奠定基礎。

Shibly等［3］對輪壤模型進行簡化，并對土壤內(nèi)聚力和內(nèi)摩擦角進行了估計。Kang［4］為解決深空探測星球土壤力學參數(shù)無法直接測量的難題，設計基于車輪沉陷量的土壤參數(shù)估計方法，基于視覺估計車輪沉陷量并進行處理，采集車輪信息，利用應用統(tǒng)計方法對月壤參數(shù)進行反演。Xue等［5］基于車輪信息設計采用最小二乘法的月壤參數(shù)識別模型，能對在軌月球車月壤力學參數(shù)進行估計，該方法需要對數(shù)據(jù)進行預處理。崔平遠等［6］通過Guass-Legendre 數(shù)值積分和Newton迭代對輪壤模型進行簡化，在線估計月壤參數(shù)。Ding等［7］建立輪壤解耦模型，通過最小二乘法預測月壤參數(shù)，但只在低滑轉率時有效。

本文設計月壤參數(shù)估計器，對月壤力學參數(shù)的識別首先對輪壤相互作用模型進行簡化，基于自適應無跡卡爾曼濾波算法對月壤參數(shù)進行識別。設置不同工況對設計的估計器效果進行了驗證。

1 月球車輪壤模型

Bernstein［8］以車輪沉陷作為出發(fā)點，得出車輪壓力和土壤沉陷量之間的關系。Bekker［9］在Bernstein研究的基礎上總結出Bekker 應力模型。Janosi［10］與Bekker正應力模型為后續(xù)研究工作奠定了理論基礎。Reece等［11］對Bekker正應力模型進行了修正，增加了最大應力角、進入角和離去角等參數(shù)，對正應力分布公式進行了完善。這些應力分布理論為后續(xù)研究人員對星球車輪壤作用力學的研究提供了基礎，并且還在不斷完善修正。

在松軟月面，車輪行駛過程中會使月壤壓縮變形，根據(jù)Bekker半經(jīng)驗法思想，輪壤之間的相互作用導致月壤發(fā)生變形可分為相互獨立的法向變形和切向變形，也即車輪對月壤的承壓特性和剪切特性，如圖1所示。

圖1 輪壤應力分布

根據(jù)Bekker 提出的正應力模型，車輪-月壤作用面的平均正壓力

式中：kc為內(nèi)聚力變形模量；b為壓板半徑；kφ為摩擦變形模量；z為壓板沉陷量；n為沉陷指數(shù)。

根據(jù)車輪沉陷模型，車輪在月壤上的靜態(tài)沉陷量

式中：W為車輪載荷；D為車輪直徑。

當車輪在滑轉時，沉陷量會有所增加，Bekker 認為滑轉沉陷量是關于滑移率的線性函數(shù)，因此車輪的總沉陷量

式中：hb為輪刺高度；s為車輪滑移率。

根據(jù)Janosi剪切應力模型，車輪下方的剪切應力τ與剪切變形j關系

式中：τmax＝（c+σtan φ）為最大剪切應力，其中，c為月壤內(nèi)聚力；φ為內(nèi)摩擦角；j0為切變模量。

根據(jù)車輪幾何關系，可得到沉陷量z和車輪進入角θ1之間關系

式中，R為車輪半徑。剪切變形

式中，θ為輪壤接觸區(qū)域角度。

車輪在實際行駛過程中，會發(fā)生滑轉，導致車輪的正壓力分布峰值并不在其正下方，而會產(chǎn)生角度的偏移

式中，c1和c2為最大應力角系數(shù)。

根據(jù)Wong-Reece的修正正應力分布模型，車輪正下方的正應力分布和剪切應力分布可分別為：

式中，θ2為車輪離去角。

車輪的垂直載荷FN、掛鉤牽引力FDP和月壤力矩TR分別為：

由式（10）～（12）構成的載人月球車輪壤相互作用模型以復雜的積分表示，很難用于求解月壤參數(shù)識別，要對輪壤作用應力-應變分布模型進行簡化，將輪地作用力進行線性化［3，12］。

令：

式中，E＝1 -exp｛-R［（θ1-θm）-（1 -s）（sin θ1-sin θm）］／j0｝。則式（8）正應力分布和式（9）剪切應力分布可分別簡化為：

對積分式求解并進行簡化，可得車輪垂向載荷FN、月壤力矩TR和掛鉤牽引力FDP表達式分別為：

式中，B＝。

由表達式可得，若要通過輪壤相互作用關系求解車輪力，還需要計算出各車輪的滑移率。

車輪滑移率

式中：v為各輪縱向輪心速度；ω為車輪轉速。

2 月壤參數(shù)估計器

不同月壤參數(shù)的變化對載人月球車牽引力、驅動力矩和垂向載荷的影響不盡相同，要確定需要識別參數(shù)。描述月壤力學特性的參數(shù)包括月壤內(nèi)聚力c、內(nèi)摩擦角φ、內(nèi)聚力變形模量kc、摩擦變形模量kφ以及切變模量j0。其中，內(nèi)聚力c對于月球車牽引和力矩的影響很小，而內(nèi)聚力c對牽引力和力矩的取值很敏感，本文不考慮對內(nèi)聚力c的識別，在參數(shù)識別過程中給定一個合理值。其他月壤力學參數(shù)內(nèi)摩擦角φ、內(nèi)聚力變形模量kc、摩擦變形模量kφ以及切變模量j0，取值變化對于車輪牽引力和力矩的影響比較明顯，對內(nèi)聚力變形模量kc和摩擦變形模量kφ統(tǒng)一考慮，令kp＝kc／b+kφ。選擇內(nèi)摩擦角φ、系數(shù)kp以及切變模量j0作為待識別參數(shù)［6］。

無跡卡爾曼濾波算法（Unscented Kalman Filter，UKF）是由無跡變換（Unscented transform，UT）和線性卡爾曼濾波算法共同構成的。無跡變換應用了高斯分布原理，其本質(zhì)是：在原始的狀態(tài)量數(shù)據(jù)中根據(jù)特定的標準來獲取部分采樣點，采樣點的均值、協(xié)方差矩陣和原始狀態(tài)數(shù)據(jù)保持一致；將其代入非線性狀態(tài)空間，可獲得映射的非線性狀態(tài)空間點集，根據(jù)原始點的選取規(guī)則，獲得UT 后的均值與協(xié)方差矩陣。根據(jù)高斯分布的原理，可獲得非線性系統(tǒng)精確到三階的UT 后的后驗均值及協(xié)方差矩陣［13］。

UKF算法步驟：

步驟1濾波初始化

式中：ˉx0為初始狀態(tài)；P0為初始協(xié)方差矩陣；E為求均值。

步驟2UT。UT 是在原始自變量數(shù)據(jù)附近按特定規(guī)則選取sigma點，使得sigma點的數(shù)據(jù)均值和協(xié)方差矩陣與原始自變量一致，即通過選取的sigma 點近似原始自變量；將采集到的sigma 采樣點數(shù)據(jù)代入非線性函數(shù)，得到因變量的均值和協(xié)方差矩陣。對于非線性變換y＝f（x），狀態(tài)量為n維，已知其均值ˉx 和協(xié)方差矩陣P，可采用UT變換來得到2n+1 個sigma點和對應的權值計算因變量的統(tǒng)計特征。

2n+1 個sigma采樣點

采樣點相應的權值

式中：n為狀態(tài)量維數(shù)；λ為縮放比例因子；α為較小正值，其決定了Sigma 點的分布狀態(tài)；β 為非負的權系數(shù)，可合并方程中高階項的動態(tài)誤差。

步驟3 根據(jù)k時刻2n+1 個sigma 采樣點計算k+1 時刻均值和協(xié)方差矩陣

式中：f為系統(tǒng)過程方程；qk為過程噪聲；Qk為過程噪聲協(xié)方差。

步驟4 將sigma 采樣點代入測量方程h得到預測觀測值

步驟5 根據(jù)預測觀測值加權求和，得到系統(tǒng)預測均值、因變量協(xié)方差以及自變量與因變量的協(xié)方差：

式中：rk為測量噪聲；Rk為測量噪聲協(xié)方差。步驟6 更新卡爾曼濾波增益

步驟7 更新系統(tǒng)狀態(tài)和協(xié)方差

盡管UKF對非線性系統(tǒng)有較高的估計精度，但存在一個問題，就是需要提前知曉系統(tǒng)噪聲統(tǒng)計特性。如果系統(tǒng)噪聲或者噪聲統(tǒng)計特性存在偏差，則可能導致系統(tǒng)估計方差偏大，使得濾波發(fā)散。

Sage-Husa 自適應無跡濾波算法（Adaptive Unscented Kalman Filter，AUKF）是在UKF 的基礎上，在遞推濾波過程中，利用采集到的信息，在時變噪聲估值器作用下，實時估計和更新非線性系統(tǒng)的過程噪聲與測量噪聲的均值和協(xié)方差，以達到減小系統(tǒng)估計誤差、增加濾波準確度和減少濾波發(fā)散［14］。使用AUKF對月球車行駛過程中系統(tǒng)過程噪聲和測量噪聲的估算步驟：

步驟1計算系統(tǒng)過程噪聲均值估計值

步驟2計算系統(tǒng)過程噪聲協(xié)方差估計值

式中，ek＝y(tǒng)k-為殘差。

步驟3計算系統(tǒng)測量噪聲均值估計值

步驟4計算系統(tǒng)測量噪聲協(xié)方差估計值

基于AUKF 對月壤參數(shù)進行參數(shù)識別，需搭建月壤參數(shù)識別系統(tǒng)的過程和測量方程。

載人月球車在月面實際行進途中，由于設定的2個估計時刻之間的間隔很短，可認為2 個估計時刻之間月壤參數(shù)維持同一個值不變，故系統(tǒng)過程方程

車輪垂向力、牽引力和力矩作為測量過程求解月壤參數(shù)，3 種力可由月球車配置的力／力矩傳感器獲得，測量方程：

3 Adams／Simulink仿真

Adams是一款系統(tǒng)動力學仿真軟件，廣泛應用在各大汽車廠商和研究院所，其仿真分析功能十分強大，涵蓋范圍廣，能滿足用戶建模的不同需求，其中Car模塊包含了一系列的車輛仿真專用模塊，可用來快速建立樣車模型，并對其各種功能指標進行評價。在Adams 建立的模型還可與其他軟件如Matlab、Solidworks、Pro／E 等進行交互，實現(xiàn)聯(lián)合仿真。在Adams建立的月球車模型如圖2 所示。

圖2 月球車Adams模型

為驗證基于AUKF月壤參數(shù)估計算法的準確性，建立Adams／Simulink聯(lián)合仿真模型，本節(jié)選擇月球車以18 km／h的恒定縱向車速分別在月壤1 和月壤2 行駛，其中估計器的3 個目標參數(shù)發(fā)生變化，其余參數(shù)不變（見表1），通過設計的AUKF 估計器可得到月壤的參數(shù)估計值，以左前輪為例，月壤參數(shù)估計器和實際值仿真結果對比如圖3 所示。

表1 兩種類型月壤參數(shù)［15］

圖3 月壤參數(shù)估計結果

由仿真結果可見，估計器設置的初值和月壤實際參數(shù)有所不同，無論是月壤1 還是月壤2，基于AUKF的月壤參數(shù)估計算法都可在0.5 s左右波動后收斂到真實值附近，且隨后一直保持穩(wěn)定。因為月壤參數(shù)估計器采用的輪壤模型對實際模型有所簡化和月球車仿真中輪壤作用關系存在一定差值，因此月壤參數(shù)估計器的估計結果和實際月壤參數(shù)之間會有一定誤差，結果收斂后的差值很小，誤差在2%以內(nèi)。

4 結語

（1）本文提出基于AUKF 的月壤參數(shù)估計方法，并進行了聯(lián)合仿真。

（2）建立月球車行駛過程中輪壤模型并進行簡化，設計AUKF 估計器，通過車輪力／力矩對月壤力學參數(shù)內(nèi)摩擦角φ、系數(shù)kp以及切變模量j0進行估計。

（3）ADAMS／Simulink 聯(lián)合仿真結果表明所設計的估計算法能夠有效識別月壤參數(shù)。