徐 明

(北京航空航天大學(xué) 宇航學(xué)院，北京100191)

譚 田 李志武

(航天東方紅衛(wèi)星有限公司 研發(fā)中心，北京 100094)

徐世杰

(北京航空航天大學(xué) 宇航學(xué)院，北京100191)

隨著人類空間活動的發(fā)展，航天器對軌道機(jī)動能力的需求越來越強(qiáng).以軌道攔截、遠(yuǎn)程交會等飛行任務(wù)為例，往往對機(jī)動初始、終止位置以及轉(zhuǎn)移時間有嚴(yán)格約束，即Lambert問題.由于Lam－bert轉(zhuǎn)移不存在解析解［1－2］，眾多學(xué)者著重利用優(yōu)化算法進(jìn)行求解:參考文獻(xiàn)［3］首先利用遺傳算法得到圓軌道之間的最優(yōu)轉(zhuǎn)移;參考文獻(xiàn)［4］通過改進(jìn)遺傳算法進(jìn)一步得到Lambert雙脈沖轉(zhuǎn)移的最優(yōu)解.但這些研究成果均建立在各種測量誤差已知的假設(shè)下;考慮到誤差分布的隨機(jī)性，所得到的最優(yōu)解并不具有普適性.

事實(shí)上，在轉(zhuǎn)移期間施加中途修正，可有效地減少交會的落點(diǎn)偏差:這在地月轉(zhuǎn)移和拉格朗日點(diǎn)轉(zhuǎn)移中已有成功應(yīng)用［5］.中途修正策略的設(shè)計(jì)包括修正次數(shù)、每次軌控時刻以及修正量;本文基于參考軌跡提出線性和非線性3種策略，并應(yīng)用Monte－Carlo和遺傳算法的聯(lián)合仿真，得到實(shí)現(xiàn)代價函數(shù)(落點(diǎn)誤差最小)在概率意義下的最優(yōu)解.

與以往直接求解Lambert問題不同，本文將限制性三體問題中求解周期性特解的微分修正算法［6－8］加以改造，得到 J2－Lambert轉(zhuǎn)移軌道并作為中途修正的參考軌跡.顯然，J2－Lambert轉(zhuǎn)移軌道基本消除攝動項(xiàng)的影響，則中途修正僅需要補(bǔ)償導(dǎo)航誤差、初始偏差修正的控制偏差等，無需補(bǔ)償攝動力等影響.

1 微分修正算法生成J2-Lambert轉(zhuǎn)移軌道

根據(jù)r0，rf及Δθ可計(jì)算二體Lambert轉(zhuǎn)移軌道［1］.顯然，考慮J2攝動后，該變軌策略將造成較大的偏差.本文基于生成Halo軌道［2］和J2不變軌道［3］的微分修正算法，設(shè)計(jì) J2－Lambert轉(zhuǎn)移軌道.

Lambert變軌策略給出位置交會的轉(zhuǎn)移軌道，控制量為位置0處的脈沖速度增量Δθ;微分修正算法將改進(jìn)ΔV及轉(zhuǎn)移時間Δt=tf－t0，以實(shí)現(xiàn)位置f處的交會.Lambert軌道轉(zhuǎn)移見圖1.

圖1 Lambert軌道轉(zhuǎn)移

貼近真值的迭代初值，可以保證微分修正算法迭代過程的收斂.由于J2攝動的數(shù)量級僅為10－3m/s2，故本文 ΔV的初值取自二體 Lambert變軌策略;為了將探測器導(dǎo)引到目標(biāo)位置，每次迭代過程都將軌道積分至Rf處，且與Rf具有相同的x(或y，z)坐標(biāo)分量，而tf即取為該軌道積分時間.顯然，0處軌道速度V0的變化ΔV將導(dǎo)致軌道積分時間tf的變化Δtf.

考查第m次迭代，對轉(zhuǎn)移末端f處的位置矢量作一階Taylor展開，可得

速度修正量ΔV應(yīng)使得

整理式(1)和式(2)，可得

實(shí)際上，微分修正算法可推廣到更高階引力場模型;但考慮到測量和軌控等誤差對落點(diǎn)偏差的影響大于非球型攝動高階項(xiàng)，即J3－/J4－Lambert轉(zhuǎn)移軌道并不能明顯改善中途修正的效果;故本文僅利用微分修正算法消除主攝動(J2)項(xiàng)引起的落點(diǎn)偏差.

2 轉(zhuǎn)移軌道的線性化

其中，單值矩陣定義為

可通過求解矩陣微分方程得到:

單值矩陣具有如下性質(zhì):

3 中途修正策略設(shè)計(jì)

3.1 基于參考軌跡的線性修正策略I

設(shè)Lambert轉(zhuǎn)移過程中進(jìn)行k次修正.根據(jù)修正次數(shù)，整個名義軌道被劃分為k+1段，如圖2所示.

圖2 Lambert轉(zhuǎn)移中的中途修正

速度增量ΔVm由于修正軌道，以實(shí)現(xiàn)探測器在時刻到達(dá)目標(biāo)位置故由式(10)可解得

該變種可以看作有k個獨(dú)立的修正策略Ⅰ組成，且每個策略僅進(jìn)行一次修正.

3.2 基于參考軌跡的線性修正策略Ⅱ

為了解決以上問題，本文擬采用隨機(jī)最優(yōu)控制理論設(shè)計(jì)修正策略Ⅱ.

為了方便表示，將δX記為Y，將 F(tm+1，tm)簡記為F(m+1，m).則探測器運(yùn)動狀態(tài)的一階近似由以下線性差分方程表示:

式中，D(m)=F(m+1，m)·［03×3I3×3］T;白噪聲向量W∈N(0，G)和N∈N(0，Q)分別來源于系統(tǒng)未建模誤差及每次修正引入的誤差和測量誤差.

交會(或攔截)位置精度與燃料消耗的加權(quán)關(guān)系，由二次型代價函數(shù)F0確定，其定義為

對于該最優(yōu)問題式(12)～式(13)，可以按照不完全狀態(tài)信息情形下離散隨機(jī)最優(yōu)控制的動態(tài)規(guī)劃法［5，9］求解.求解過程分為以下 3 步:

稱為Kalman增益.

線性最優(yōu)濾波器:

從而最優(yōu)速度修正為

3.3 策略Ⅰ和Ⅱ的關(guān)系

不考慮燃料代價，則代價函數(shù)F0中加權(quán)矩陣取 GR=I3×3，K→∞，即 K－1=0.通過求解 Riccati方程及Kalman增益，可得

式(23)表明，策略Ⅰ本質(zhì)上是策略Ⅱ的特殊形式.

3.4 基于Lambert迭代的非線性修正策略Ⅲ

微分修正即可設(shè)計(jì)參考軌跡，還可設(shè)計(jì)修正策略.整個Lambert轉(zhuǎn)移過程依次在…，k時刻進(jìn)行k次修正;對于任意修正點(diǎn)，探測器的狀態(tài)記為和，完成該轉(zhuǎn)移的剩余時間為則根據(jù)第1節(jié)發(fā)展的“J2－Lambert轉(zhuǎn)移軌道的微分修正生成算法”，可以解算出轉(zhuǎn)移到期望終端狀態(tài)Rf和Vf的修正脈沖ΔVm.

顯然，策略Ⅲ獨(dú)立于名義軌道，具有更強(qiáng)的魯棒性;策略Ⅰ(特別是其變種)，僅是策略Ⅲ的一階近似，故策略Ⅲ具有更高精度.策略Ⅰ和Ⅱ僅需進(jìn)行線性運(yùn)算，且反饋矩陣經(jīng)離線設(shè)計(jì)并裝訂在星上處理器以減少計(jì)算量;相反地，策略Ⅲ需要進(jìn)行非線性的迭代運(yùn)算，這將耗費(fèi)相當(dāng)?shù)挠?jì)算資源.

4 最優(yōu)修正時刻選擇

完整的修正策略應(yīng)包括修正量輸出和修正時刻確定.對于給定修正量輸出策略，希望設(shè)計(jì)合適的修正點(diǎn)，以實(shí)現(xiàn)代價函數(shù)在概率意義下的最優(yōu).

而修正時刻的確定，需要針對具有不同統(tǒng)計(jì)特性的誤差散布來源，對整個隨機(jī)過程建立復(fù)雜的統(tǒng)計(jì)學(xué)模型;且中途修正不同于連續(xù)控制:對于整個動力系統(tǒng)的演化過程，僅需有限次的不連續(xù)控制.因此，很難給出概率意義下最優(yōu)修正時刻的解析形式，本文將基于Monte－Carlo和遺傳算法的聯(lián)合仿真給出數(shù)值求解方法.

記εi為最后落點(diǎn)誤差，其中i=90% ～99%，表示該誤差的置信水平.則對于任意的修正時刻τ，存在映射 Θ:τ→εi(τ).顯然，映射 Θ 存在于統(tǒng)計(jì)學(xué)意義，并且可由Monte－Carlo仿真得到確定.

本文期望找到合適的τ，以使得εi達(dá)到最小值.因此，最優(yōu)修正時刻的選擇轉(zhuǎn)化為映射Θ的最優(yōu)化計(jì)算;而工程上，置信水平一般取為90%，95%，99%，即該問題本質(zhì)上是多目標(biāo)優(yōu)化問題.幸運(yùn)地是，遺傳算法的計(jì)算結(jié)果表明:該多目標(biāo)最優(yōu)化模型的Pareto最優(yōu)解為同一結(jié)果，即存在多目標(biāo)優(yōu)化最優(yōu)解.

為了優(yōu)化算法的實(shí)現(xiàn)，需對優(yōu)化變量修正時刻τ作必要的離散化處理:對整個轉(zhuǎn)移時間Δt=tf－ t0等間隔離散為 N 段，即└t0，t1，…，tN－1，tf┘;并離線計(jì)算出各子區(qū)間的單值矩陣F(tm+1，tm)，m=0，1，…，N －1，用于構(gòu)造修正策略Ⅰ或Ⅱ.

5 仿真算例

根據(jù)表1的軌道根數(shù)確定Lambert轉(zhuǎn)移軌道，如圖3所示;整個轉(zhuǎn)移時間Δt=9 000 s.在J2項(xiàng)攝動的影響下，該轉(zhuǎn)移軌道將有65.09 km的落點(diǎn)偏差;在21×21階引力攝動下，落點(diǎn)偏差將達(dá)到65.134 km;經(jīng)過“生成 J2－Lambert轉(zhuǎn)移軌道的微分修正算法”矯正后，該轉(zhuǎn)移軌道在21×21階引力攝動下偏差僅為35.9243 m.

表1 初始軌道和目標(biāo)瞬時軌道根數(shù)

由于存在導(dǎo)航誤差E1和轉(zhuǎn)移軌道的初始偏差E2以及Lambert脈沖速度偏差E3，將會造成探測器真實(shí)飛行軌道偏離名義軌道，故中途修正是必要的;而中途修正將會引入修正速度偏差E4.

以下誤差源認(rèn)為服從Gauss分布，均值都為零，方差各為:E1和E2的位置偏差200 m，速度偏差0.01m/s(1σ);E3和E4的速度大小偏差0.1%，速度在徑向方向偏差0.573°(3σ).而速度法平面內(nèi)的偏差滿足［0，2π］上的均勻分布.

本文采用修正策略Ⅰ來制導(dǎo)Lambert轉(zhuǎn)移，并利用“Monte－Carlo和遺傳算法的聯(lián)合仿真”求解最優(yōu)修正時刻.遺傳算法的參數(shù)取為:種群數(shù)為40，遺傳代數(shù)為25，選擇代溝為0.9，交叉概率為0.7，變異概率為 0.0017;Monte－Carlo 仿真次數(shù)為500.優(yōu)化變量τ的離散間隔N取為300.

優(yōu)化目標(biāo)為尋找ε90%，ε95%和ε99%的最小值，而遺傳算法的數(shù)值結(jié)果顯示:該多目標(biāo)最優(yōu)化模型的Pareto最優(yōu)解為同一結(jié)果，即存在多目標(biāo)優(yōu)化最優(yōu)解 τ*=390 s，對應(yīng)的落點(diǎn)誤差分別為8.82 km(置信度90%)，9.94 km(置信度95%)和10.82 km(置信度99%).最優(yōu)解的尋優(yōu)過程見圖4，而最優(yōu)修正策略的落點(diǎn)誤差分布見圖5.

圖3 Lambert轉(zhuǎn)移軌道

圖4 遺傳算法的尋優(yōu)過程

圖5 最優(yōu)修正策略的落點(diǎn)誤差分布

6 結(jié)論

本文應(yīng)用Monte－Carlo法和遺傳算法的聯(lián)合仿真求解Lambert轉(zhuǎn)移中途修正的全局概率最優(yōu)策略.以往利用優(yōu)化算法直接求解Lambert問題，需要已知各種測量誤差的大小，即所得到的最優(yōu)解并不具有普適性.為了有效地減小交會的落點(diǎn)偏差，可在轉(zhuǎn)移期間施加中途修正.本文針對轉(zhuǎn)移的最終落點(diǎn)誤差，設(shè)計(jì)3類中途修正策略以適應(yīng)不同的精度需要.以第Ⅰ類修正策略為例，應(yīng)用Monte－Carlo和遺傳算法的聯(lián)合仿真，得到在概率意義下的最優(yōu)修正時刻.

中途修正的參考軌跡設(shè)計(jì)，基于限制性三體問題中求解周期性特解的微分修正算法構(gòu)造，以消除攝動項(xiàng)的影響，則中途修正僅需要針對導(dǎo)航誤差、初始偏差修正的控制偏差等進(jìn)行補(bǔ)償.

本文應(yīng)通過Monte－Carlo法和遺傳算法的聯(lián)合仿真發(fā)現(xiàn)“不同置信水平具有相同的概率最優(yōu)策略”結(jié)論;該結(jié)論對于設(shè)計(jì)工程實(shí)用的修正策略具有重要意義，而如何從理論上證明上述結(jié)論將是下一步的工作重點(diǎn).

References)

［1］Battin R H.An introduction to the mathematics and methods of astrodynamics［M］.New York:AIAA Inc，1987

［2］韓潮，謝華偉.空間交會中多圈Lambert變軌算法研究［J］.中國空間科學(xué)技術(shù)，2004(5):9－14 Han Chao，Xie Huawei.Research on algorithm of loopy Lambert transfer in space rendezvous［J］.Chinese Space Science and Technology，2004(5):9 －14(in Chinese)

［3］Spencer D B，Kim Y H.Optimal spacecraft rendezvous using genetic algorithms［J］.Journal of Spacecraft and Rockets，2002，39(6):859－865

［4］陳統(tǒng)，徐世杰.基于遺傳算法的最優(yōu) Lambert雙脈沖轉(zhuǎn)移［J］.北京航空航天大學(xué)學(xué)報(bào)，2007，33(3):273 －274 Chen Tong，Xu Shijie.Optimal Lambert two－impulse transfer using genetic algorithm［J］.Journal of Beijing University of Aeronautics and Astronautics，2007，33(3):273 －274(in Chinese)

［5］Xu Ming，Xu Shijie.Trajectory and correction maneuver during the transfer from Earth to Halo orbit［J］.Chinese Journal of Aeronautics，2008，21(3):200 －206

［6］Xu Ming，Xu Shijie.Nonlinear dynamical analysis for displaced orbits above a planet［J］.Celestial Mechanics and Dynamical Astronomy，2008，102(4):327 －353

［7］Xu Ming，Xu Shijie.Structure－preserving stabilization for Hamiltonian system and its applications in solar sail［J］.Journal of Guidance，Control and Dynamics，2009，32(3):997 －1004

［8］Xu Ming，Xu Shijie.J2invariant relative orbits via differential correction algorithm［J］.Acta Mechanica Sinica，2007，23(5):585－595

［9］Noton M.Spacecraft navigation and guidance［M］.London:Springer，1998