龍 也，劉一武

(1．北京控制工程研究所，北京100190;2．空間智能控制技術(shù)重點實驗室，北京100190)

0 引言

以實現(xiàn)0．1 km定點著陸精度為目標的火星精確著陸問題是目前火星探測進入、下降和著陸(Entry，Descent and Landing，EDL)過程制導(dǎo)、導(dǎo)航與控制技術(shù)研究的熱點［1－2］。而作為EDL過程中氣動環(huán)境最惡劣、高度跨度最大，參數(shù)變化及不確定性最多的階段［1］，火星進入段的末端狀態(tài)精度對前述目標實現(xiàn)非常關(guān)鍵。為此需深入研究影響進入段制導(dǎo)精度的具體因素及高精度魯棒制導(dǎo)的實現(xiàn)途徑［2］。

研究表明，進入狀態(tài)偏差是引起進入段末端狀態(tài)偏差的主要擾動來源之一［3］。目前，在“火星科學(xué)實驗室”等實際探測任務(wù)中得到應(yīng)用的“解析預(yù)測校正進入制導(dǎo)方式”存在一定缺陷，它對進入狀態(tài)偏差敏感，需嚴格保證精度，否則初始偏差經(jīng)進入段誤差傳播可能引起較大著陸誤差［2］，而相關(guān)魯棒制導(dǎo)律的研究［4－5］則集中在制導(dǎo)跟蹤的魯棒性上，缺乏針對擾動來源的直接優(yōu)化。相對而言，基于脫敏最優(yōu)控制(Desensitized Optimal Control，DOC)的制導(dǎo)方法［6－9］提供了一種新的設(shè)計思路，它通過在性能指標中添加特定形式的“系統(tǒng)對狀態(tài)擾動敏感度”罰項來獲得具有相應(yīng)脫敏特性的標稱軌跡，從而實現(xiàn)針對性的制導(dǎo)魯棒性提升。但目前該類方法通常采用反轉(zhuǎn)邏輯［10－11］控制橫程偏差，以傾側(cè)角反轉(zhuǎn)可瞬時完成為假設(shè)前提，未考慮實際過程中傾側(cè)角調(diào)整能力的限制。

本文在考慮前述問題的基礎(chǔ)上針對“進入狀態(tài)偏差”進行三自由度脫敏軌跡設(shè)計研究。全文安排如下:第一節(jié)介紹動力學(xué)模型、現(xiàn)有縱向脫敏方法及非脫敏方法，第二節(jié)分析傾側(cè)角調(diào)整性能受束即非理想反轉(zhuǎn)條件下反轉(zhuǎn)邏輯對制導(dǎo)精度的影響，說明縱向脫敏存在的局限性并指出對于脫敏類制導(dǎo)方法采用三自由度設(shè)計的必要性，第三節(jié)討論三自由度模型下敏感度傳播奇異問題并給出本文設(shè)計，第四節(jié)進行仿真分析。

1 預(yù)備知識

1．1動力學(xué)模型與基本約束

考慮自轉(zhuǎn)影響的三自由度火星進入段動力學(xué)方程如下:

式中:r為探測器到火星中心距離，v為速度，γ為航跡傾角，s為航程，θ為經(jīng)度，λ為緯度，ψ為航跡方位角，ωmrs為火星自轉(zhuǎn)角速度，g=μ/r2為探測器所處位置的重力加速度，μ為火星引力常數(shù)，σ為傾側(cè)角。L和D分別為氣動升力和阻力加速度:

式中:ρ為火星大氣密度，采用標準指數(shù)型大氣模型ρ=ρsexp［－(h/hs)］，h=r－rM為所處位置高度，rM為火星半徑，ρs為火星表面大氣密度，hs為火星大氣尺度高度，CL和CD分別為氣動升力和阻力系數(shù)，Sr為探測器參考面積，m為探測器質(zhì)量。此外也可采用如下方法計算升阻力加速度:

式中:B為彈道系數(shù)，Γ為升阻比系數(shù)。

受限于探測器物理參數(shù)極限及開傘安全要求，進入段需考慮的基本約束包括最大過載約束、最大動壓約束、最大熱流率約束等三項過程約束及末端開傘動壓約束、速度約束等兩項末端約束:

式中:nmax，qmax，，qfmax，vfmax等為前述約束的最大容許值，kq，rn，N，M為熱流率計算相關(guān)系數(shù)。

1．2縱向脫敏方法

縱向脫敏方法采用縱向簡化動力學(xué)方程，由式(1)～式(4)構(gòu)成并忽略自轉(zhuǎn)，記為如下形式:

式中:x=［r，v，γ，s］T為系統(tǒng)狀態(tài)向量，u=cosσ為控制輸入，t為時間，f(x，u)為式(1)～式(4)右側(cè)項式。t∈［t0，tf］時，令x(t)=X(t|t0，x0)表示其在初始條件x(t0)=x0下的解，則敏感度矩陣S(t|t0，x0)按下式定義［9］:

簡記為S(t|t0)，它滿足如下敏感度傳播方程［7］:

跟蹤控制律:

式中:u*(t)和x*(t)分別為t時刻的標稱控制量和標稱系統(tǒng)狀態(tài)，K(t)為狀態(tài)反饋增益陣，η(t)為反饋增益乘子，通常η(t)=1，但考慮控制飽和問題時，部分文獻采用如下形式［8］:

按式(19)積分就能得到軌跡各處對于進入狀態(tài)偏差的敏感度陣，進而可根據(jù)不同設(shè)計目標［6，9］計算敏感度罰項，通過優(yōu)化求解過程得到具有相應(yīng)脫敏特性的標稱軌跡。

1．3非脫敏方法

非脫敏方法用于設(shè)計簡單的標稱軌跡，本文中作為參照衡量其它兩種方法性能，軌跡同樣通過優(yōu)化求解過程得到，動力學(xué)方程及約束與縱向脫敏方法一致，目標函數(shù)為:

式中:J1=－h(huán)(tf)，Jn2=u2，cn0為懲罰因子。

2 非理想反轉(zhuǎn)對制導(dǎo)精度的影響

反轉(zhuǎn)邏輯的優(yōu)點在于偏差控制策略的簡便性及橫縱程設(shè)計的解耦特性，因此得到廣泛應(yīng)用［9－13］。但需注意的是，實際過程中探測器的傾側(cè)角調(diào)整能力有限，特別是弱機動能力航天器的調(diào)整能力更差，這使得真實反轉(zhuǎn)過程不能滿足假設(shè)前提，從而可能對制導(dǎo)精度造成影響。目前部分文獻［12－13］已考慮到該類問題，但通常僅關(guān)注如何補償偏差，較少對其影響制導(dǎo)精度的具體途徑進行分析。

事實上，非理想反轉(zhuǎn)引起的直接變化主要包括以下兩個方面:一是影響垂直升力分量，使得控制指令偏離標稱值，二是導(dǎo)致制導(dǎo)過程出現(xiàn)開環(huán)飛行狀態(tài)。本節(jié)從這兩點出發(fā)討論其對制導(dǎo)精度的影響，并對縱向脫敏方法在非理想反轉(zhuǎn)條件下的真實受擾程度進行分析(圖表均采用第4節(jié)設(shè)計的相應(yīng)標稱軌跡測試得到)。

2．1升力分量偏離標稱值的影響

由傾側(cè)角反轉(zhuǎn)基本定義可知反轉(zhuǎn)過程中探測器垂直升力分量相對標稱值偏大，因此航跡傾角會不可避免的有所增加。單次反轉(zhuǎn)過程額外增加的γ總量可按下式估計［14］:

其中單次反轉(zhuǎn)時長Δt可通過假定傾側(cè)角按恒定的最大變化角速度反轉(zhuǎn)進行估算:

式中:σ0為反轉(zhuǎn)前σ瞬時值。同時垂直升力增量

令LT=CLqSr，將式(25)代入式(23)中整理得:

2．2反轉(zhuǎn)過程制導(dǎo)開環(huán)的影響

反轉(zhuǎn)過程中制導(dǎo)邏輯默認不進行狀態(tài)閉環(huán)［12］，而弱機動能力航天器反轉(zhuǎn)耗時較長，這使得其制導(dǎo)過程會較長時間處于開環(huán)狀態(tài)，導(dǎo)致偏差消除能力下降，間接削弱了制導(dǎo)魯棒性。

由圖2和表1可知，縱向設(shè)計的“非脫敏軌跡”應(yīng)用于真實制導(dǎo)過程時反轉(zhuǎn)次數(shù)多達八次，反轉(zhuǎn)過程耗時74．73 s，約占總時長30%左右。由圖3可知，制導(dǎo)末端高度及速度標準差會隨傾側(cè)角調(diào)整能力的下降而趨于增大，即系統(tǒng)對進入狀態(tài)偏差的敏感度上升，制導(dǎo)魯棒性下降。

圖1 非理想反轉(zhuǎn)條件下開環(huán)標稱制導(dǎo)末端狀態(tài)偏差與的關(guān)系Fig．1 Relation between open-loop nominal guidance terminal state error and under non-ideal inverse condition

2．3縱向脫敏受擾程度分析

圖3 非理想反轉(zhuǎn)條件下閉環(huán)制導(dǎo)末端狀態(tài)標準差與的關(guān)系Fig．3 Relation between close-loop guidance terminal state standard deviation andnder non-ideal inverse condition

為充分反映非理想反轉(zhuǎn)對縱向脫敏方法實際制導(dǎo)過程的影響，本節(jié)采用“縱向脫敏軌跡”和“非脫敏軌跡”進行對照，觀察它們在無進入狀態(tài)偏差條件下的閉環(huán)制導(dǎo)過程以及考慮進入狀態(tài)偏差問題時的末端狀態(tài)偏差抑制能力。由圖4和圖5對比可知，采用縱向脫敏方法時系統(tǒng)在反轉(zhuǎn)完成后易因狀態(tài)偏差過大而出現(xiàn)持續(xù)的控制飽和。由圖6和圖7對比可知，非理想反轉(zhuǎn)條件下“非脫敏軌跡”末端高度偏差范圍為±1 km，速度偏差范圍為±50 m/s，而“縱向脫敏軌跡”末端高度偏差范圍為 －8～2 km，速度偏差范圍為－250～50 m/s。

由此可見，對于縱向脫敏方法，非理想反轉(zhuǎn)容易引起持續(xù)的控制飽和并會導(dǎo)致其所設(shè)計軌跡在蒙特卡洛仿真測試中不但不能達到預(yù)期的制導(dǎo)魯棒性增強效果，反而出現(xiàn)末端狀態(tài)偏差散布程度數(shù)倍于非脫敏方法的嚴重失效問題。

圖4 非脫敏軌跡制導(dǎo)過程控制剖面對照Fig．4 Control profile comparison of non-desensitized trajectory guidance process

綜上所述，非理想反轉(zhuǎn)對弱機動能力航天器制導(dǎo)精度影響明顯，且對縱向脫敏方法影響尤為嚴重，這說明采用反轉(zhuǎn)邏輯進行橫向偏差控制的傳統(tǒng)制導(dǎo)策略并不適用于弱機動能力航天器精確制導(dǎo)。因此對“脫敏類制導(dǎo)方法”來說，取消反轉(zhuǎn)邏輯，考慮三自由度模型下的脫敏設(shè)計十分必要。

圖5 縱向脫敏軌跡制導(dǎo)過程控制剖面對照Fig．5 Control profile comparison of longitudinal desensitized trajectory guidance process

圖6 非理想反轉(zhuǎn)條件下非脫敏軌跡末端狀態(tài)偏差分布Fig．6 Terminal state error distribution of non-desensitized trajectory under non-ideal inverse condition

圖7 非理想反轉(zhuǎn)條件下縱向脫敏軌跡末端狀態(tài)偏差分布Fig．7 Terminal state error distribution of longitudinal desensitized trajectory under non-ideal inverse condition

3 三自由度脫敏設(shè)計

動力學(xué)方程包括式(1)～式(7)并考慮傾側(cè)角調(diào)整動態(tài):

式中:σv為傾側(cè)角角速度，uc為傾側(cè)角角加速度。優(yōu)化求解過程以uc為系統(tǒng)輸入。

3．1控制及狀態(tài)約束

基礎(chǔ)約束包括式(12)～式(16)，并附加傾側(cè)角范圍約束、角速度及角加速度約束:

附加末端狀態(tài)約束:

末端航程stf為自由狀態(tài)，經(jīng)緯度與預(yù)期值θtf，λtf一致時標稱軌跡結(jié)束。

3．2敏感度傳播方程

首先，類似于縱向脫敏方法，考慮x為完整的狀態(tài)量［r，v，γ，s，θ，λ，ψ］T，f(x，u)為式(1)～ 式(7)右項。此時，敏感度傳播方程中?f/?u項為:

顯然u=1時?f/?u最后一項變?yōu)闊o窮，出現(xiàn)“奇異問題”，這會導(dǎo)致傳播方程無法繼續(xù)推演，因此前述x與f(x，u)的選取并不合理。

應(yīng)注意到，在僅考慮縱向剖面的現(xiàn)有脫敏方法中x=［r，v，γ，s］T包括了對末端開傘安全最為重要的高度、速度項，而航程s相當于將耦合的橫縱程脫敏問題轉(zhuǎn)化為單一量表示。因此，敏感度傳播方程仍采用式(19)，跟蹤控制律中η(t)=1－(u*(t))2。

3．3目標函數(shù)

目標函數(shù)的選取包括末端開傘高度和軌跡脫敏兩項。一方面，開傘高度是保證探測器安全著陸的重要因素，且開傘高度越高意味著探測任務(wù)點可選范圍更大［1］。另一方面，S(tf|t0)包含末端狀態(tài)對進入狀態(tài)偏差的敏感度信息。因此目標函數(shù)J設(shè)為:

式中:J1=－h(huán)(tf)為開傘高度罰項，c0≥0為敏感度懲罰因子，J2為進入狀態(tài)偏差敏感度項:

式中:ci≥0，i=1，2，3，4為權(quán)重因子，表示對各狀態(tài)脫敏特性的注重程度。

4 仿真分析

本節(jié)首先給出仿真參數(shù)，其次分別采用“非脫敏方法”，“縱向脫敏方法”與“三自由度脫敏方法”設(shè)計標稱軌跡，最后通過蒙特卡洛仿真校驗三自由度脫敏的改進效果。其中部分軌跡設(shè)計及仿真結(jié)果已在前文引用，不再重復(fù)給出。

4．1仿真參數(shù)

反轉(zhuǎn)邏輯采用的航向誤差走廊參數(shù)如下［11］:

式中:stgo為待飛航程，航向誤差Δaz超出最大允許航向誤差Δazm時傾側(cè)角進行反轉(zhuǎn)。其它設(shè)計及仿真參數(shù)見表2～表5。

表2 相關(guān)物理參數(shù)Table 2 Related physical parameters

此外 σt0=0，=0，K=－［0．01，0．005，50，0．001］。

對于前兩種設(shè)計方法，表4中前四項約束保留，且stf=735 km，smax=stf并有umin=－0．8，umax=0．8。

4．2軌跡設(shè)計結(jié)果

軌跡的優(yōu)化求解均采用高斯偽譜法［15］實現(xiàn)。圖8～圖9為三自由度脫敏軌跡設(shè)計結(jié)果。

表3 控制及過程約束參數(shù)［1］Table 3 Control and process constraint parameters［1］

表4 邊界條件及約束參數(shù)Table 4 Boundary conditions and constraint parameters

表5 蒙特卡洛仿真參數(shù)Table 5 Monte-Carlo simulation parameters

圖8 c0=0．02、h(tf)=10210 m時三自由度脫敏軌跡標稱傾側(cè)角剖面Fig．8 Nominal bank angle profile of three degree of freedom desensitized trajectory corresponding to c0=0．02 and h(tf)=10210 m

圖9中傾側(cè)角加速度與速率軌跡表明本文方法在設(shè)計時已保證所設(shè)計軌跡在探測器傾側(cè)角調(diào)整能力范圍之內(nèi)。

圖9 三自由度脫敏軌跡完整信息Fig．9 Complete information of three degree of freedom desensitized trajectory

4．3蒙特卡洛仿真校驗

圖10為“縱向脫敏軌跡”在縱向剖面內(nèi)校驗時的蒙特卡洛仿真結(jié)果，圖11為“三自由度脫敏軌跡”在三自由度模型下校驗時的蒙特卡洛仿真結(jié)果。

圖10表明“縱向脫敏軌跡”的末端高度偏差為±50 m，速度偏差為±15 m/s。圖11表明“三自由度脫敏軌跡”的末端高度偏差為±10 m，速度偏差為 ±1 m/s，航程偏差為 ±0．2 km，縱向偏差為 ±0．5 km，橫向偏差為±4 km。由此可知，本文設(shè)計相比現(xiàn)有縱向脫敏方法可更大程度降低系統(tǒng)對進入狀態(tài)偏差的敏感程度，具有顯著增強的魯棒性能。

圖10 在縱向剖面內(nèi)校驗時縱向脫敏軌跡的蒙特卡洛測試結(jié)果Fig．10 Monte-Carlo test results of longitudinal desensitized trajectory when tested in longitudinal profile

圖11 三自由度脫敏軌跡的蒙特卡洛測試結(jié)果Fig．11 Monte-Carlo test results of three degree of freedom desensitized trajectory

5 結(jié)束語

現(xiàn)有縱向脫敏方法通常采用反轉(zhuǎn)邏輯控制橫程偏差，但本文研究表明非理想反轉(zhuǎn)會引起末端狀態(tài)偏差、使系統(tǒng)對進入狀態(tài)偏差敏感度上升，并且導(dǎo)致弱機動能力航天器采用縱向脫敏方法制導(dǎo)時出現(xiàn)末端狀態(tài)偏差散布程度數(shù)倍于非脫敏方法的失效情況。因此，對于火星進入段制導(dǎo)，本文放棄采用反轉(zhuǎn)邏輯，針對進入狀態(tài)偏差問題提出了三自由度脫敏設(shè)計。仿真結(jié)果表明，該方法有效解決非理想反轉(zhuǎn)條件下縱向脫敏方法的嚴重失效問題，且末端狀態(tài)偏差抑制能力更強。

目前本文僅考慮了進入狀態(tài)偏差問題，針對“氣動參數(shù)不確定”、“過程氣動力擾動”的脫敏設(shè)計及同時考慮三類擾動形式的完整脫敏方案有待進一步研究。

［1］Braun R D，Manning R M．Mars exploration entry，descent and landing challenges［J］．Journal of Spacecraft and Rockets，2007，44(2):310－323．

［2］ 崔平遠，胡海靜，朱圣英．火星精確著陸制導(dǎo)問題分析與展望［J］．宇航學(xué)報，2014，35(3):245－253．［Cui Ping-yuan，Hu Hai-jing，Zhu Sheng-ying．Analysis and prospect of guidance aspects for Mars precision landing［J］．Journal of Astronautics，2014，35(3):245－253．］

［3］ 李爽，彭玉明，陸宇平．火星EDL導(dǎo)航、制導(dǎo)與控制技術(shù)綜述與展望［J］．宇航學(xué)報，2010，31(3):621－627．［Li shuang，Peng Yu-ming，Lu Yu-ping．Review and prospect of Mars EDL navigation guidance and control technologies［J］．Journal of Astronautics，2010，31(3):621－627．］

［4］Mooij E，Mease K D，Benito J．Robust re-entry guidance and control system design and analysis［C］．AIAA Guidance，Navigation，and Control Conference and Exhibit，Hilton Head，USA，August，2007．

［5］Lu W M，Bayard D S．Guidance and control for Mars atmospheric entry:adaptivity and robustness［R］．AIAA Paper，1999．

［6］Shen H J，Hans S，Richard W P．Desensitizing the pin-point landing trajectory on Mars［C］．AIAA/AAS Astrodynamics Specialist Conference and Exhibit，Honolulu，Hawaii，August 18－21，2008．

［7］Seywald H，Kumar R R．Desensitized optimal trajectories［J］．Advances in the Astronautical Sciences，1996，93(1):103－115．

［8］Seywald H．Desensitized optimal trajectories with control constraints［J］．Advances in the Astronautical Sciences，2003，114:737－743．

［9］Li S，Peng Y M．Mars entry trajectory optimization using DOC and DCNLP［J］．Advances in Space Research，2011，47(3):440－452．

［10］Harpold JC，Graves Jr CA．Shuttle entry guidance［J］．Journal of the Astronautical Sciences，1979，27(3):239－268．

［11］Tu K Y，Munir M S，Mease K D，et al．Drag-based predictive tracking guidance for Mars precision landing［J］．Journal of Guidance，Control，and Dynamics，2000，23(4):620－628．

［12］Carman G，Ives D，Geller D．Apollo-derived Mars precision lander guidance［C］．AIAA Atmospheric Flight Mechanics Conference and Exhibit，Boston，MA，August 10－12，1998．

［13］Brunner C W．Skip entry trajectory planning and guidance［D］．Iowa:Iowa State University，2008．

［14］Mooij E，Space D．Model reference adaptive guidance for reentry trajectory tracking［C］．AIAA Guidance，Navigation，and Control Conference and Exhibit，Rhode Island，USA，August 16－19，2004．

［15］Garg D，Michael P，Hager W，et al．A unified framework for the numerical solution of optimal control problems using pseudospectral methods［J］．Automatica，2010，46(11):1843－1851．